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Abstract 



The achievement of sustained nuclear fusion in magnetically confined plasma relies 
on efficient confinement of alpha particles, which are high-energy ions produced by 
the ftision reaction. Such particles can excite instabilities in the frequency range of 
Alfven Eigenmodes (AEs), which significantly degrade their confinement and threat- 
ens the vacuum vessel of future reactors. In order to develop diagnostics and control 
schemes, a better understanding of linear and nonlinear features of resonant inter- 
actions between plasma waves and high-energy particles, which is the aim of this 
thesis, is required. In the case of an isolated single resonance, the description of 
AE destabilization by high-energy ions is homothetic to the so-called Berk-Breizman 
(BB) problem, which is an extension of the classic bump-on-tail electrostatic problem, 
including external damping to a thermal plasma, and collisions. A semi-Lagrangian 
simulation code, COBBLES, is developed to solve the initial-value BB problem in 
both perturbative (Sf) and self-consistent (full-/) approaches. Two collision models 
are considered, namely a Krook model, and a model that includes dynamical fric- 
tion (drag) and velocity-space diffusion. The nonlinear behavior of instabilities in 
experimentally-relevant conditions is categorized into steady-state, periodic, chaotic, 
and frequency-sweeping (chirping) regimes, depending on external damping rate and 
collision frequency. The chaotic regime is shown to extend into a linearly stable region, 
and a mechanism that solves the paradox formed by the existence of such subcriti- 
cal instabilities is proposed. Analytic and semi-empirical laws for nonlinear chirping 
characteristics, such as sweeping-rate, lifetime, and asymmetry, are developed and 
validated. Long-time simulations demonstrate the existence of a quasi-periodic chirp- 
ing regime. Although the existence of such regime stands for both collision models, 
drag and diffusion are essential to reproduce the alternation between major chirp- 
ing events and quiescent phases, which is observed in experiments. Based on these 
findings, a new method for analyzing fundamental kinetic plasma parameters, such 
as linear drive and external damping rate, is developed. The method, which consists 
of fitting procedures between COBBLES simulations and quasi-periodic chirping AE 
experiments, does not require any internal diagnostics. This approach is applied to 
Toroidicity-induced AEs (TAEs) on JT-60 Upgrade and Mega- Amp Spherical Toka- 
mak devices, which yields estimations of local kinetic parameters and suggests the 
existence of TAEs relatively far from marginal stability. The results are validated by 
recovering measured growth and decay of perturbation amplitude, and by estimating 
collision frequencies from experimental equilibrium data. 



Resume 



Le succes de la fusion nucleaire par confinement magnetique repose sur un confine- 
ment efficace des particules alpha, qui sont des ions hautement encrgctiqucs produits 
par les reactions de fusion. De telles particules pcuvent exciter des instabilites dans 
le domaine de frequence des modes d'Alfven (AEs) qui degradent leur confinement 
et risquent d'endommagcr I'enceinte a vide de reacteurs futurs. Afin de developper 
des diagnostiques et moyens de controle, une meilleure coniprchcnsion des comporte- 
ments lineaire et non-lineaire des interactions resonantes entre ondes plasma et par- 
ticules energetiques, qui constitue le but de cette these, est requise. Dans le cas 
d'une resonance unique et isolee, la description de la destabilisation des AEs par des 
ions energetiques est homothetique au probleme de Berk-Breizman (BB), qui est une 
extension du probleme classique de I'instabilite faisceau, incluant un amortissement 
externe vers un plasma thermiquc. ct des collisions. Un code semi-Lagrangien, COB- 
BLES, est developpe pour resoudre le probleme aux valeurs initiales de BB selon 
deux approches, perturbative (Sf) et auto-coherente (full-/). Deux modeles de col- 
lisions sont consideres, a savoir un modele de Krook, et un modele qui indue la 
friction dynamique et la difi^usion dans I'espace des vitesses. Le comportement non- 
lineaire de ces instabilites dans des conditions correspondantes aux experiences est 
categorise en regimes stable, periodique, chaotique, et de balayage en frequence (sif- 
fiet), selon le taux d'amortissement externe et la frequence de collision. On montre que 
le regime chaotique deborde dans une region lineairement stable, et Ton propose un 
mecanisme qui resoud le paradoxe que constitue I'existence de telles instabilites sous- 
critiques. On developpe et valide des lois analytiques et semi-empiriques regissant les 
caracteristiques non-lineaires de siflilet, telles que la vitesse de balayage, la duree de vie, 
et I'asymetrie. Des simulations de longue duree demontrent I'existence d'un regime de 
sifflets quasi-periodiques. Bien que ce regime existe quel que soit I'un des deux modeles 
de collision, la friction et la diffusion sont essentielles pour reproduire I'alternance 
entre sifflets et pcriodcs do rcpos. telle qu'observee experimentalcmcnt. Grace a 
ces decouvertes, on developpe une nouvelle methode pour analyser des parametres 
cinetiques fondamentaux du plasma, tels que le taux de croissance lineaire et le taux 
d'amortissement externe. Cette methode, qui consistc a faire correspondrc les sim- 
ulations de COBBLES avec des txpcriences d'AEs qui presentent des sifflets quasi- 
periodiques, ne requiert aucun diagnostique interne. Cette approche est appliquee a 
des AEs induits par la toroidicitc (TAEs) sur les machines JT-60 Upgrade ct Mega- 
Amp Spherical Tokamak. On obtient des estimations de parametres cinetiques locaux 
qui suggerent I'existence de TAEs relativement loin de la stabilite marginale. Les 
resultats sont valides en recouvrant la croissance et dccroissance de I'amplitudc des 
perturbations mesurees, et en estimant les frequences de collision a partir des donnees 
experiment ales d'equilibre. 
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Chapter 1 

Introduction 



This manuscript has been written in an effort toward pedagogy, privileging clarity 
over details, physical pictures over mathematical developments. We hope it can be 
used as an introduction to the Berk-Breizman problem and its numerical computa- 
tion. However we cannot pretend being able to introduce notions such as nuclear 
fusion, magnetic plasma confinement, Magneto-HydroDynamics (MHD), kinetic the- 
ory, tokamak geometry, or particle orbits, with as much clarity as in well-known 
textbooks. For this reason, the reader is assumed to be familiar with the basis of the 
above fields. In this introduction we expose the background and motivation for our 
work. 

1.1 Energetic particles-driven AEs 

In an ignited tokamak operating with a deuterium-tritium mix, the confinement of a- 
particles is critical to prevent damages on the first- wall and to achieve break-even. The 
reason is that these high energy particles must be kept long enough in the plasma 
core to allow enough of their energy to heat thermal populations by slowing-down 
processes. A major concern is that high energy ions can excite plasma instabilities in 
the frequency range of Alfven Eigenmodes (AEs), which significantly enhance their 
transport. Ever since the recognition of this issue in the 1970's, considerable progress 
has been made in the theoretical understanding of the principal Alfvenic instabilities. 
However, the estimation of the mode growth rate is complex, and the question of 
their stability in ITER [ACH+Ol] remains to be clarified. In addition, estimation of 
the effect on transport and development of diagnostics and control schemes require 
significant progress on our understanding of nonlinear behaviors. 

We limit our framework to the tokamak configuration. In this work, we focus our 
interest on the Toroidicity induced Alfven Eigenmode (TAE) |CCC85) . which is a 
shear-Alfven wave (SAW) located in a gap of the SAW continuum that is created by 
toroidal coupling of two successive poloidal modes, and which can be destabilized by 
resonant interactions with high-energy ions. TAEs have been observed to be driven by 
g-partic les in burnin g plamas [WSB+QSllNEB+OTj . lon- Cyclotron-Resonance Heating 
(ICRH) |WWC+94| . and Neutral Beam Injection (NBI) [WFP+91llWSD+91) . In this 
work we consider only the latter situation (NBI-driven TAEs). 

1.2 The BB problem as a paradigm for TAEs 

In general, TAEs are described in a three-dimensional (3D) configuration space. How- 
ever, near a phase-space surface where the resonance condition is satisfied (resonant 
phase-space surface), it is possible to obtain a new set of variables in which the 
perturbation is described by a one-dimensional (ID) effective Hamiltonian in 2 con- 
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ugated variables |BBP97al lBBP+97cl IWB98| IGDPN+OS] . if we assume an isolated 



J 

single mode. In this sense, the problem of kinetic interactions between TAEs and 
fast-particles is homothetic to a simple ID single mode bump-on-tail instability. The 
so-called Berk-Breizman (BB) problem |BBY93[ lBBP97al IBBP+ 97cl IBBP96| is a gen- 
eralization of the bump-on-tail problem, where we take into account an external wave 
damping accounting for background dissipative mechanisms, and a collision operator. 
Observed qualitative and quantitative similarities between BB nonlinear theory and 
both global TAE simulations |WB98[ |PBG+04| and experiments tFBB+98[ IHFSOO] 
are an indication of the validity of this reduction of dimensionality. Similar arguments 
exist for other Alfven wave instabilities such as the geodesic acoustic mode (GAM) 
IBBB+OS] . and the beta Alfven eigenmode (BAE) |Ngu09| . These analogies enable 
more understanding of fully nonlinear problems in complex geometries by using a 
model that is analytically and numerically tractable. This approach is complemen- 
tary to heavier 3D analysis. 



1.3 Frequency sweeping 

A feature of the nonlinear evolution of AEs, the frequency sweeping (chirping) of 
the resonant frequency by 10-30% on a timescale much faster than the equilib- 
rium evolution, has been observed in the plasma core region of tokamaks JT-60 
Upgrade (JT-6 0U) |KKK+99] . DIII-D 'Hei95], the SmaU Tight Aspect Ratio Toka- 



mak (START) jMGS+99l . the Mega Amp Spherical Tokamak (MAST) PBG+Oi) , 
the National Spherical Torus Experiment (NSTX) |FGB06| . and in stellerators such 
as the Large Helical Device (LHD) jOYT+06] . and the Compact Helical Stellerator 
(CHS) |TTT+99] . In general, two branches coexist, with their frequency sweeping 
downwardly (down-chirping) for one, upwardly (up-chirping) for the other. In many 
experiments, asymmetric chirping has been observed, with the amplitude of down- 
chirping branches significantly dominating up-chirping ones. Chirping TAEs have 
been reproduced in 3D simulations with a hybrid MHD/drift-kinetic code |TSTI03j . 
and with a drift-kinetic perturbative code |PBG"'"04] . Qualitatively similar chirping 
modes are spontaneously generated by the BB model, and have been shown to corre- 
spond to the evolution of holes and clumps in the velocity distribution. Theory relates 
the time evolution of the frequency shift with the linear drive 7^ and the external 
damping rate 7^ BBP97bj , when the evolution of holes and clumps is adiabatic. The 
idea of using chirping velocity (or sweeping rate) as a diagnostic for growth rates 
looks promising. If we assume that the mode is close to marginal stability, 7^ ~ 7^, 
and if we further assume that holes and clumps are in the adiabatic regime, then 
the growth rates are simply obtained by fitting analytic prediction to experimental 
measurement of chirping velocity. However, these two assumptions are not verified in 
general, which limits the applicability of this simple approach. In most experiments, 
chirping events are quasi-periodic, with a period in the order of the millisecond. On 
the one hand, the long-time simulation of repetitive chirping with 3D codes is a heavy 
task, which has not been undertaken yet. On the other hand, simulations of the BB 
model with many chirping events have been performed [VBSCOTl ILBSIO) , but such 
solutions do not feature any quasi-periodicity. In this sense, repetitive chirping that 
qualitatively agrees with an experiment has not been reproduced, neither with ID 
nor 3D simulations. 



1.4 Aim of this thesis 

Our main goals are to improve our understanding of wave-particle nonlinear resonant 
interactions, develop diagnostics and identify control parameters for AEs. From these 
backgrounds, a straight-forward plan is to: 
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• clarify the link between BB model and AEs, 

• extend BB theory where it is insufficient to explain experimental observations, 

• analyze experimental data by applying the BB model to AEs. 

1.5 Outline 

In brief, we reduce the problem of TAEs to a numerically and analytically tractable 
paradigm, the BB problem. We make a survey of linear and nonlinear physics of the 
BB model, and work on improving our understanding, by extending theory and by 
using systematic numerical observations, focusing on the frequency sweeping regime. 
Armed with new findings and robust numerical tools, we finally go back to the original 
problem of TAEs, explaining quantitative features of experimental observations, and 
developing a new method for accurate linear predictions. 

In Chap. [2j we review the physics of the TAE, and simplify the description of 
the problem from 6D to 2D in phase-space, around a single resonant phase-space 
surface. The first step in this procedure is to isolate the gyromotion, which is de- 
coupled from the wave for typical TAEs with / ~ 100 kHz, by changing variables 
to the guiding-center coordinates. This change of variables is based on the so-called 
Lie transform perturbation theory, which we review in order to get a better grasp of 
involved hypothesis. Then we show how to reduce the guiding-center Hamiltonian, 
as well as collision operator and background mechanisms, and discuss limitations of 
this simplified description. In Chap. [3] we review basic nonlinear physics of the BB 
problem. We develop and verify numerical tools that we use in following chapters. 
In particular, we develop and verify a kinetic code to solve the initial-value problem. 
We show that numerical simulations in experimentally relevant conditions, with cold- 
bulk and weak, warm-beam distributions, require a careful choice of advection scheme. 
Among the family of so-called Constrained Interpolation Profile (CIP) schemes, a lo- 
cally conservative version (CIP-CSL) displays the best convergence properties. As an 
intermediate summary, we cast the analogies between BB model and ID description 
of TAEs. In Chap. |4j we investigate kinetic features of self-consistent interactions 
between an energetic particle beam and a weakly unstable electrostatic wave in ID 
plasma, within the BB framework, with an emphasis on chirping regime and sub- 
critical instabilities. We show that the nonlinear time-evolution of chirping in ID 
simulations can be used to retrieve information about linear input parameters with 
good precision. We identify a regime where chirping events are quasi-periodical. This 
regime exists whether the collision model is annihilation/creation type, or takes into 
account dynamical friction and velocity-space diffusion. Based on these findings, in 
Chap. [5] we propose a new method to estimate local linear drive, external damping, 
and collision frequencies from the spectrogram of magnetic field variations measured 
by a Mirnov coil at the edge of the plasma. This method, which relies on a fitting 
of normalized chirping characteristics between the experiment and BB simulations, 
is described and applied to TAE experiments in MAST and JT-60U. We show that 
the BB model can successfully reproduce features observed in the experiment only 
if the collision operator includes drag and diffusion terms. We find a quantitative 
agreement for the diffusion frequency, and a qualitative agreement for the drag fre- 
quency, between the values obtained with our fitting procedure, and an independent 
estimation obtained from experimental equilibrium measurements. In the conclusion 
(Chap, [g]), we summarize our findings, and we suggest numerical and experimental 
investigations these findings enable. 
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Chapter 2 

Energetic particle-induced 
TAEs 



The TAE is a shear- Alfven wave located in a gap of the SAW continuum that is created 
by toroidal coupling of two successive poloidal modes, which can be destabilized by 
energetic particles. This chapter deals with the modelization at several levels of com- 
plexity of an isolated single-mode TAE. In Sec. |2.1[ we give a short review of the basic 
TAE physics, and of the state-of-the-art of its linear stability analysis, emphasizing a 
need for improved accuracy. When the TAE is an isolated single-mode, it is possible 
to reduce the problem to a simple harmonic oscillator, by expanding the perturbed 
Hamiltonian around a resonant phase-space surface. This reduction from 3D to ID 
becomes particularly transparent in angle- action variables (0;^, Ji), once the pertur- 



bation has been put in the form exp(iZiQ;i — iwi), where U are integers. In Sec. 2.2 
we show how to obtain the latter form. The first step is to separate the gyromotion, 
which does not resonate with typical TAEs, by changing variables to guiding-center 
coordinates. This procedure can be done while conserving the Hamiltonian properties 
by making use of Lie transform perturbation theory, which we review thoroughly. In 
Sec. |2.3[ we show how to reduce the Hamiltonian, collision operator, and background 
damping mechanisms from their 3D description to a ID model. 

2.1 Review 

This short review is intended to introduce notions and notations that we use when 
applying the BB model to the TAE, and to motivate our linear analysis of TAEs. For 
a more comprehensive review of energetic-particle driven AEs, see Ref. |Hei08| . 

2.1.1 Toroidal geometry 

Hereafter, we make use of magnetic flux coordinates [RDH03] , also called as straight 
field-line coordinates, (r, 0, ^), where r, 9 and C are minor radius, poloidal and toroidal 
angle, respectively, to describe the nested poloidal flux surfaces of the equilibrium 
magnetic field Bo- In these coordinates. Bo takes the following form, 

Bo = Vx X VM-C), (2.1) 

where x(r) is the poloidal fiux (divided by 27r), and the safety factor, defined by 

is the flux surface-averaged number of toroidal rotation that a field line undergoes in 
one poloidal rotation. 
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Fig. 2.1: Alfven continuum for n = 1, with and without coupUng between m — 2 and 
m — Z poloidal modes. The q profile has been modeled by q(r) = 1.2 + 2.1 (r ja)^-^ ^ 



and the electronic density by n^(r) = 0.11 + 1.57(1 — r /a ) 
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Other 



parameters are Bq = 1.2 T, Rq = 3.3 m, and a = 0.96 m. Note that the discrepancy, 
relatively far from resonance, between the upper branch of coupled continuum and 
the uncoupled m — 2 branch, is accounted by errors of second-order in the aspect 
ratio. Similar discrepancy is observed in Fig. 1 of Ref. |FD89) for example. 



2.1.2 Physics of the TAE 

In an homogeneous magnetized plasma, linear ideal-MHD arguments show the exis- 
tence of a shear- Alfven wave of frequency wa with the dispersion relation 



2 



where 



1 2 2 
fc|| Wa, 



So 



(2.3) 



(2.4) 



is the Alfven velocity, and fcy is the wave number in the direction of the equilibrium 
magnetic field Bo- Let us consider axisymmetric toroidal plasmas. In the cylindrical 
limit, the periodicities of the system require that there exists two integers, a toroidal 
mode number n and a poloidal mode number m, such that 



fell = 



n — m/q{r) 
Rf) 



(2.5) 



where i?o is the distance from the symmetry axis of the tokamak to the magnetic 
axis. In a non-homogeneous plasma in a sheared magnetic field, both fcii and ua are 



functions of r. The simple dispersion relation Eq. (2.3) is still valid in this configura- 
tion ICC86j , and it is called the Alfven continuum. Since phase velocity is a function 
of radius, a wave packet with finite radial extent would suffer from phase-mixing, the 
so-called continuum damping. Except for energetic particle modes, which are outside 
of the scope of this work, resonant drive by fast particles is not enough to overcome 
this damping. However, a toroidal coupling of two successive poloidal modes m and 
m+1 breaks up the continuous spectrum. This is illustrated in Fig. |2.H which shows 
the Alfven continuum for n = 1, m = 2 and m -I- 1 = 3, in cylindrical geometry, 
where two poloidal continuum are decoupled, and in toroidal geometry, with a two- 
mode coupling model |CC86[ IFD89) . The latter is obtained with equilibrium plasma 
parameters corresponding to JT-60U shot E32359 at i = 4.2s, which we analyze in 
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Chap. [5j assuming concentric circular magnetic flux surfaces, retaining toroidicity 
efi^ects in the first order in inverse aspect ratio. Though we show only the w > 
half-plane, the continuum is symmetric with respect to u. Coupled modes are (n, m) 
and (— n, — m — 1) for > 0, and (n, m + 1) and (— n, —to) for w < 0. The gap 
is centered at a radius 7'a such that q{rx) = (to + l/2)/n, where the two continu- 
ous spectra would cross in the absence of coupling, and where |A;||| — l/2qRQ. The 
resulting discrete eigenmode is a TAE, at a frequency uja = vx/2qRQ. 

For a deuterium plasma with typical magnetic field Bq ^ IT and density n.i ~ 
102°TO-3, the Alfvenic energy is Ea = miv\/2 ^ lOfceF, which is in the range of 
passing particles induced by NBI. For ITER parameters, Ea ~ IMeV, which is in 
the range of passing a-particles born from fusion reactions. In both cases, TAEs can 
be driven unstable by resonance with energetic particles. For far-passing particles, 
the resonance condition is f2 = uja, where 

il = nuj(^ + liL)g, (2-6) 

where — W||/i?o find ujg = v\\/qRo are frequencies of toroidal motion and poloidal 
motion, respectively, and I = —to for co-passing particles, I = m for counter-passing 
particles |TS98j . Since we analyze TAEs driven by co-injected ions, we can simplify 
following discussions by considering only co-passing particles. Then, the resonance 
condition is 

Vn Vn 

WA - n — + TO—- = 0. (2.7) 
ito qKo 



2.1.3 Linear stability 

In theory, the linear growth rate 7 ~ 7l — is positive when the drive by fast 
particles overcomes damping processes to background plasma. The growth rate can 
be estimated either by linear stability co des, such a s PENN |JAVV95j . TASK/WM 
\FKm\ . NOVA-K |Che92) , or CASTOR-K [BBB+02| ; or by gyr o- or drift- k inetic per - 
turbative nonlinear initial value codes, such as FAC |CBB+97| or HAGIS |PAC+98| . 
The analysis requires internal diagnostics that are not always available. 

The linear drive depends on several factors such as spatial and energy gradients 
of EP distribution, and alignment between particle orbit and the eigenmode. In 
the limit of large aspect ratio, analytic t heory |FD89| yields successful estimations 
of 7l for well-defined numerical models [TSW"'"95] . However, in the general case, 
complicated factors cited above forbid accurate analytic estimation, as yet. Moreover, 
improvements in measurement of EP distributions are needed to allow estimation of 
7i for experimental TAEs. 

The external damping rate 7^ involves complicated mechanisms, which include 
continuum damping (since parts of the mode extend spatially into the continuum) 
|ZC92j . radiative damping |MM92| . Landau damping with thermal species [BF921 
IZCS96] . and colhsional damping by trapped electrons |GS92| . Most of these mech- 
anisms are still under investigation, and cannot be quantified by existing theory. 
Experimentally, 7^ can be estimated by active measurements of externally injected 
perturbations |FBB+95l IFBB+OO] . However, the applicability of this technique is 
limited to dedicated experiments. 

Moreover, if the system is close to marginal stability, 7 is sensitive to small varia- 
tions of driving and damping terms, and is very sensitive to plasma parameters such 
as the q profile. In addition, the existence of unstable AEs in a regime where linear 
theory predicts 7 < 0, or subcritical AEs, has not been ruled out. To access the 
subcritical regime, nonlinear analysis is necessary. Therefore, accurate linear stabil- 
ity analysis requires significant theoretical and experimental improvements, or a new 
approach. 
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2.2 Angle-action description 



2.2.1 Kinetic description 

Since the kinetic equation is at the center of interest of this thesis, it is useful to get 
a deep understanding of its meaning, by reviewing its derivation from fundamental 
principles. The steps we summarize here are detailed in Ref. |KT86| for example. 

A many-body system is completely described by the microscopic distribution 
N{z,t) — ~ ^i)' where Zi = {xi,Vi), xi and vi are the position and ve- 

locity of a particle indexed as i, and the sum is taken over all particles. To simplify 
the present discussion, we consider a single-species system of Np particles, without 
external forces, and normalize the total phase-space volume, assumed constant, to 1. 
Substituting Newton equations of motion into the partial time derivative of N yields 
the so-called Klimontovich-Dupree equation, 

dN dN ■ dN 

-+..-+a— .-=0. (2.8, 

where a™^'^^°{z) is the acceleration due to microscopic electromagnetic fields, at the 
exception of those due to a particle located at z, if such a particle exists. The 
microscopic electromagnetic fields obey Maxwell equations, where density and current 
are velocity moments of N . 

Since it is impossible to reproduce any many-body experiment at the microscopic 
level, it is much more efficient to take an ensemble point-of-view, where distributions 
and fields are smooth functions of phase-space. The statistical properties are com- 
pletely determined by the distribution Fn{z'^,z'2^ . . . , z'j^ , t), where F/vdVi . . . AVn^ 
is the probability of finding, at a time t, particle 1 within dVi, particle 2 within dV2, 
. . . and particle Np within dV/Vp , and dV^ is an infinitesimal phase-space volume at 
the neighborhood oi z = z'^. By integrating F^v over all z[ but z^, we can define the 
one-particle distribution function /i, where fi{z'^,t)AVi is the probability of finding 
one particle within dVi at t. fi is related to the microscopic distribution by 

/,,..,«, = <M. (2.8) 

where (N) is the ensemble average J FpfNdz[ . . . dz'^ . Similarly, by integrating 
F/v over all z'^ but z{ and z'2, we can define the two-body distribution function 
f2{z[, z'2,t), which is related to the microscopic distribution by 

. f) - iN{zi,t)N{z2,t)) 5{z^~Z2) 

j2(Zi,Z2,t) - — fi(Zi,t), (2. 10) 

Np 

where we used a large particle number approximation, Np ^1. In the absence of 



atomic and nuclear processes, Np is a constant, then the ensemble average of Eq. (2.8) 
yields 

dh ^ ^.a/i ^ ^.a/i ^ d/i 

dt dx dv dt 



(2.11) 



coll. 



where the average acceleration a = ^a™"^""") is given by electromagnetic fields that 
obey Maxwell equations, where density and current are velocity moments of /i. The 



collision term, 

d/i 
dt 



is shown to vanish in the absence of particle interactions, yielding Vlasov equation. 



which is valid on a time-scale much shorter than a coUisional time-scale. Eq. (2.11) 
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Fig. 2.2: Illustration of Liouville's theorem. Time-evolution of an infinitesimal volume 
in a 2D phase-space {x, v), delimited by a solid curve at t = and by a dashed curve 
at i = t2. The phase-space volume is constant, dVi = dV2, and, in the absence of 
collision, the number of particles inside the volume is constant too. 



is not a closed equation, because the second part of the collision term involves ex- 
pressions of the form (NN). Under the Coulomb approximation, which forbids any 
retardation effect, and which is valid if the thermal velocity is much slower than the 
speed of light, we can reduce the latter term as a function of /2. Similarly, the equa- 
tion which gives the evolution of /2 involves terms of the form (NNN), and so on. 
Altogether, we have a chain of equations for which we need a closure. A collision 
operator is an approximative statistical operator that accounts for particle interac- 
tions, which provides such closure. A clear review of collision operators is given in 
Ref. |HS02| . In this thesis, we focus on a Fokker-Planck collision operator, which 
is based on the fact that, in a Tokamak plasma, collisions are dominated by binary 
Coulomb collisions with small-angle deflection. 

In the following, we write /i simply as /. The kinetic equation, Eq. (2.11), can 
be put in Hamiltonian form. 



(2.13) 



coll. 



where h is the Hamiltonian, and {} are Poisson brackets. The l.h.s. is the variation 
of / following particle orbits, or so-called Lagrangian derivative of /, noted D(/. 

The fact that, in the absence of collisions, / is conserved along particle orbits, can 
be seen as a direct consequence of Liouville's theorem, which states that the density in 
phase-space is constant along particle orbits. Let us consider an infinitesimal volume 



of phase-space d Vi surrounding a particle at zi at t = ti, as illustrated in Fig. |2.2 
The particle and all points of dVi evolve following the equations of motion until a 
time t = t2 where the particle is at Z2, and dVi is changed to a volume dV2. Liouville's 
theorem arises from the following properties, 



Since the points of the boundary of dV follow the equations of motion, in the 
absence of collisions, the number of particles inside dV is constant; 



• Time-evolution can be seen as a series of infinitesimal canonical transformations 
generated by the Hamiltonian; 

• Poincare's invariant: the volume element in phase-space is a canonical invariant. 
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2.2.2 Review of Lie transform theory 

When working in canonical variables, it is easy to exploit the properties of a Hamilto- 
nian system. However, for some systems it is difhcult to find canonical variables, and 
the most convenient variables may not be canonical. Moreover, if some variables are 
canonical for the unperturbed system, they may be non-canonical for the perturbed 
one. By applying the Lie near-identity transformation theory to the phase-space La- 
grangian, one can change any Lagrangian into a simpler form in coordinates that 
reveal the symmetries of the system, and use this formulation to study the properties 
of a perturbed Hamiltonian system in any arbitrary noncanonical variables |CL83| . 



Noncanonical Hamiltonian Mechanics 

Consider a Hamiltonian system with N degrees of freedom. Hamilton equations are 
given by the application of a variational principle to the scalar Lagrangian L, in some 
arbitrary coordinates on the 2N + 1-dimensional extended phase space, {t, z^), 



6 J i(z^(A),A)dA = 0, (2.14) 

where A is an arbitrary parameter. Hereafter, Greek indices a, ^ and v run from to 
2N, whereas Latin indices a and b run from 1 to N, and i, j, k from 1 to 2N. Thus, 
any Hamiltonian system is characterized by its Lagrangian 

L = 7.^, (2.15) 

or equivalently by its fundamental one-form, 7^dz^. In the canonical extended phase- 
space z^ = {t,q°',p°-), when = t is the time coordinate, 

70 = -/^(^^^), (2.16) 
7a = Pa, (2.17) 

la+N = 0, (2.18) 

where h is the Hamiltonian. However, in noncanonical variables, all the 7^ may be 
nonzero. 



From Eq. (2.14), which has the same form in any new coordinate system with 



7/jd2:^ = r^dZp, Euler-Lagrange equations are obtained as 

dz" 



where w is a tensor defined by 



^^u^ = 0, (2.19) 



and its restriction to the symplectic coordinates is the Lagrange tensor. Thus 
the flow Axz^ is an eigenvector of oj^jy with eigenvalue 0. The solution is unique only 
after a normalization. In physical terms, we can take X = t. Since the Jacobian of 
the coordinate transformation 7^(2'^) = (— /i,pa,0) is nonsingular, we can invert the 
Lagrange tensor. The Poisson tensor J — then yields the equations of motion, 

( ^7i ^70 



dzo " I az" dz- ' ■ ^^-^^^ 
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Lie transform technique 



We consider a fundamental one-form 7 = 70 + £71 which consists of an equihbrium 
part 7o for which the flow is well-known, and a small perturbation 71. We want 
to study the effect of the perturbation on the flow. The strategy is to search for a 
near-identity transformation that will reveal the symmetries of the perturbed system. 
Indeed, if the fundamental one-form is independent of some coordinate z°' , then, as 
an application of Noether's theorem. 



dz" d^a dz'' dz^ 



= 0, (2.22) 



revealing an exact invariant. The general form of a near-identity transformation with 
a small parameter e is 

Z^" = z^' + eZ'^jiz) + e'^Z^^.{z) + ... (2.23) 

Rather than an expansion in e which is difficult to invert, we express the transfor- 
mation in operator form. We denote a forward transformation Z^ = Z^[z,e), and 
a backward transformation z'^ = Zj^{Z,e). In the Lie transform technique, the 
coordinate transformation is specified by a generator 5^ such that 

dZ^ 

-g^iz,e) = g^{Zf{z,e)), (2.24) 

and Zl{z,0) = zf". 

The forward transformation of a scalar /(z) into F(Z, e) is given by 

F = e-'^^f, (2.25) 
/ = e+'^oF, (2.26) 

where Lg is defined by its action on a scalar Lgf — g'^d^f , and its action on a 
one-form (Lgj)^ = g'^ {d^lii — 9^i7i/)- The transformation of coordinates is just a 
special case of scalar transformation, 

Z"^ = g-'-^f/", (2.27) 

where /"(z) = z" = Z^{Z) is the coordinate function. The transformation of a 
one-form 7(2) into T(Z,e) is given by the functional relationships, 

r = e~"'^s7 + dS, (2.28) 
7 = e+'^^T + ds, (2.29) 

where S and s are scalar functions. 



Higher order perturbation theory 



We now consider a fundamental one-form in the form of an expansion 7 = 70 + £71 + 
e^72 + ■ . .. We introduce a push-forward transformation operator T — . . .T3T2T1, 
where each T„ = e~'^ ^" is a Lie transform operator, and L„ is a short for Lg^ . 
Substituting these new definitions into Eq. (2.28), we have F = T7 + dS', which yields 
for each successive order, Tq = 70, and 



Ti = d5i - Lijo + 71 

r2 = d52 - ^270 + 72 



r2 

,■^^170 



dS„ 



for n ^ 0. 



(2.30) 
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Let us recall that our aim is to simplify the fundamental one-form in the new 
coordinates F. To this aim we have to solve successively the latter equations for the 
gauge scalar Sn and the generating vector g!^ (Each C„ is given by 7„ and the result 
of the preceding order n — 1). In these equations, the generating vector appears only 
in the one-form (L„7o)^ = g!^LOo^^. We already discussed that wo^ij/ has a unique 
null eigenvector. Then we can add any multiple of this eigenvector to g!^ without 
changing the equations we are now trying to solve. Let us suppose this eigenvector 
has a nonzero component in the time direction, then we can set 



5„ = 0, 



(2.31) 



so that the time-coordinate doesn't change {Z'^ = ^ t). At this point, we are left 
with 27V -I- 1 unknowns, namely 2N generating functions gj, and one scalar gauge 5'„. 
A priori we should be able to bring the 2 -I- 1 components of r„ into a simpler form. 
A good strategy is to make all its symplectic components Tni vanish by choosing 



9n — (diSn + Cni) Jq'- 

Let us now focus on the time component of the new one-form. 



(2.32) 



LnO 



-Hri 



dtSn 



9n^0j0 + C'nO- 



(2.33) 



It is convenient to introduce the lowest order velocity vector, defined as the time 
derivative along the unperturbed orbits of the coordinates : 



dz^ 
"dT 



(2.34) 



Substituting the unperturbed equations of motion (2.21 ), we find that the scalar gauge 



is given by its total time derivative over the unperturbed orbit. 



Dt 



(2.35) 



In integrating the latter equation, we want to avoid any secularity effect. Then we 
should remove any secular perturbation by taking 



(2.36) 



where the average is to be taken over the unperturbed orbits. 



Finally, in the new coordinates = [e+'^^f/'^] z, the new one-form is given by 

F„ = {V^Cr,^,)^ dt, (2.37) 
if we choose the generating functions 

gi = {d,Sn + Cm) (2.38) 
where the gauge scalar is given by 

Sn = - jv^nt., (2.39) 

and a tilde represents the oscillatory part. 

To illustrate the benefits of Lie transform theory compared to classical perturba- 
tion theory, and to quantify its validity limit, we apply it to a simple mathematical 
problem in App. [X] 
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2.2.3 Guiding-center Lagrangian 

Guiding-center theory provides reduced equations of motion of a charged particle in 
a slowly-varying (in time and space) electromagnetic field, where the fast gyromotion 
is decoupled from the relatively slow drifting motion of the guiding-center. Guiding- 
center theory is based on the closeness to the limit of a fixed and uniform magnetic 
field, where the orbit of a particle in a frame following its gyration center is a circle. 
The perturbation from this situation is quantified by a small parameter e, if we assume 
the following ordering, 

e - w/r^c - P/L < 1, (2.40) 

where Qc — eB/m is the cyclotronic frequency (or gyrofrequency) , p is the Larmor 
radius (or gyroradius), a; is a characteristic frequency of interest, and L is the scale- 
length of field variation. 

On the one hand, the standard derivation of guiding-center equations of motion 
is based on an averaging procedure |Nor63| . which removes important properties of 
the Hamiltonian formulation. On the other hand, the modern derivation [Lit83j is 
based on Lie-transform perturbation theory. This approach preserves the validity 
of Noether's theorem and the validity of Liouville's theorem, to each order in e. In 
addition, expansion to arbitrary order is straightforward. Moreover, since we keep 
the Hamiltonian formulation, further reductions of dimensionality are still possible 
with Lie transforms. 

The starting point is the one- form in canonical cartesian phase-space {x, p), 

7 = p-dx - h{x,p,t)dt, (2.41) 

where p — mv + eA, and A is the vector potential. The one-form can be expressed 
in noncanonical phase-space {x,v), and written as an expansion in e, 

7 = 70 + £71 + (2.42) 

where 

70 — eA{x,t) • dx — eipo{x,t), (2.43) 

71 = mv-dx - (^eipi{x,t) + ^v^^ dt, (2.44) 

and where ip is the scalar potential. 

In Ref. |Lit83] . Lie-transform theory is applied to change variables to new coor- 
dinates, where the new one-form F does not depend on the gyroangle. When this 
procedure is carried up to the second order, we obtain the guiding-center one-form, 

r = {mUb{R,t) + eA{R,t)) -dR + Md^ - Hdt, (2.45) 

where 

H = eip{R,t) + + {e/rn)MB{R,t), (2.46) 

and the guiding-center coordinates Z = {R, U, M, ^) as 



R = 


X - e—a, 


(2.47) 


U EE 




(2.48) 


M EE 




(2.49) 




tan-if"''0> 


(2.50) 


\v-e-ij 



where (ei,e2,b) is an orthogonal unit vectors system, a = cos(^)ei — sin(^)e2, and 
c = a X b. 
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In Eq. (2.45), the angle-action variables of the gyromotion, M), appear in 
canonical form in the symplectic part, and H does not depend on ^. Thus, since the 
gyromotion is irrelevant to fast-ion/TAE interactions, we can drop the term Md^ in 
the guiding-center one-form. 

2.2.4 Guiding-center action-angle variables 

In Ref. |MH90] , the zeroth order guiding-center one-form Tq is put in the form, 

To = Ped9 + P^dC + MAS, - H{Pe,PQ,e,M)At, (2.51) 

where Pq is the toroidal canonical angular momentum, 

P^ = -ex{r) + m^RaV\\, (2.52) 

and Pe is the poloidal canonical angular momentum. The latter expression can be 
used as a starting point to develop an action-angle formalism where the perturbed 
Hamiltonian takes a standard form. In Ref. |BBP95b] . a canonical transformation is 
performed to obtain 

To = JedB + JqAQ + J^d^ ~ H{Jg, J^;, J^)dt, (2.53) 

where ^ = f^c, = we, and ( — uj(^ are the unperturbed frequencies of the gyromotion, 
poloidal motion, and toroidal motion, respectively. 9 and C reduce to the geometric 
angles 6 and C if we neglect finite aspect ratio effects (But we do not neglect them, 
since we would remove the toroidicity from which the TAE originates). 



2.2.5 Application to TAE 

Although TAEs resonate mainly with passing particles, when the source of high- 
energy ions is isotropic, a large fraction of the energy transfer may be accounted by 
resonance with the bounce-motion (or banana motion) of toroidally trapped particles 
|TS98| . However, for tangential NBI ions, to which we confine our analysis, it is 
sufficient to describe resonance with far passing particles. 

In a gauge where the perturbed scalar potential is zero, the TAE can be described 
by a perturbation Hamiltonian, 

Hi = -eAi±-Vgc, (2.54) 

where v^^ is the guiding-center velocity. Here, Ai± is the perpendicular part of the 
perturbed vector potential, and we have neglected a second order term Ai± • Ai±. 
In a small plasma pressure (small /3) limit, we can neglect parallel gradients, then 
Ai_L can be split into magnetic compression, 

A5_L = bo X VH, (2.55) 

and magnetic shear, 

^i_L = V$ - (bo-V$)6o, (2.56) 

where bo = Bq/Bq. For the TAE, which is a shear Alfven wave, the latter part only 
is relevant, hence the excitation is described by a single scalar function $. 
In Ref. |BBP95b] the perturbed Hamiltonian is put in the form 

Hi = V{J) cos{l-a - tut), (2.57) 

in arbitrary tokamak geometry, where (a, J) are angle-action variables for the solvable 
unperturbed Hamiltonian Hq(J), and I — ('1,^2,^3) is a triplet of three integers. 
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Substituting Eq. (2.56) into Eq. (2.541 yields 



Hi 



-eVgc-V* + e (bo • V<I>) 60 • t'g 



d^ 



e(6o- V$)w|i 



(2.58) 
(2.59) 



where d$/dt|Q is the Lagrangian derivative of $ along unperturbed particle orbit, 
which can be removed from the Lagrangian without altering Euler equations. 

For a single n/m mode of frequency oj, the eigenfunction takes the following form, 



m+l 



$(r,6',C) = C(i)e-*"*-*^(*) ^ $,(r)e"'^-*"' + c.c. 



(2.60) 



l—m 



where C and ip are the amplitude and phase of the wave. Substituting the eigenfunc- 
tion and the expression Bq = {dx/dr)Vr X V{q9- C) into Eq. ( |2.59[ ), the TAE 
excitation is reduced to 



m+l 



Hi = -ieC{t)e 



-lUJt — lip{t) 



[uj - wybo • ("VC - mVe)] 



c.c, 



l—m 



(2.61) 

where we have neglected time-derivation of slowly varying phase and amplitude of the 
wave. ^ ^ 

Then we change the variables to the canonical angle-actions (0, Jg, Jq), in order 
to express the perturbed Hamiltonian as a Fourier series in 9, 



-t-00 



Hi = -zeC(t)e-'"'-''^We™^ ^ VpiJeJcY^^ + c.c. 



(2.62) 



where 



d6 
2^' 



m+l 



-ipO—inC, 



$i(r)e"'^-''^ [uj - W||bo- (nVC-mV6l)] 



(2.63) 



Formally, the problem is expressed in the desired form of Eq. (2.57), but the nu- 
merical computation of the Fourier components Vp, which is needed for quantitative 
comparison of absolute physical quantities between 3D and ID model, requires to 
relate geometric angles and canonical angles, which may be complicated depending 
on the equilibrium configuration. However, since each particle of the resonant phase- 
space surface interact with the wave in the same way at the same frequency (though 
with different strength), comparison of quantities that are normalized to the mode 
frequency is possible, even without evaluating Vp coefficients. 



2.3 Reduction to a one dimensional problem 

If we consider only small toroidal mode numbers n, n and n + 1 modes are isolated. 
Let us consider a single toroidal mode number. On the one hand, since on a flux 
surface r — r,„ where a poloidal mode number m is centered, the safety factor is 
Qi^m) = (2771 -|- l)/(2n), then we can estimate the distance Ar — r,„+i — r„i between 
two neighbouring m modes by writing Arq' « qi^m) ~ <l{Tm+i), as 

Ar « — . (2.64) 

nq' 

On the other hand, the characteristic width of TAE modes 5r is of the order of 
5r ^ r^/nqRo |CCC85| . Hence, for typical parameters, 6r/Ar ~ {rm/Ro)S ^ 1, 
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where 5* = rq' /q is the magnetic shear. Thus, TAEs have a two-scale radial structure, 
the larger scale corresponding to the enveloppe of the TAE. In our analysis, we assume 
that the number of m harmonics involved is small enough to consider resonances one 
by one, as isolated n, m mode. The latter hypothesis is reasonable for sufhciently 
core-localised, low-n TAEs. We must keep in mind, though, that high-n TAEs are 
likely to be destabilized in future devices such as ITER, in which case it may be 
necessary to take into account multiple-wave resonances. 

The evolution of the distribution f{x,v,t) of energetic ions is described in 3D 



configuration space by the kinetic equation (2.13), with the perturbed Hamiltonian 



Eq. (2.571. In the following, we reduce the problem to a ID Hamiltonian system, by 



considering a single n, m mode. 

2.3.1 Reduction of the Hamiltonian 

Formally, the resonant phase-space surface, J — {Jr such that Jj^^ = F{Jiii, Jr2)}, 
is defined by a function F. The resonance condition, 

l-n{jR) = w, (2.65) 

where f2(J) = jrj'{J)i is satisfied on the resonant phase-space surface. 



Once the perturbed Hamiltonian has been put in the form of Eq. (2.57), we can 
reduce the problem to one action and one angle |Lic69[ IGDPN+OS] , by performing 
a canonical transformation J ■ da — Hdt = I ■ dip — H'dt + dS with the generating 
function 

S = -I-ijj + hil-a - ujt) + hai + /aaa + F{hj2)a3- (2.66) 

This procedure yields, 

Ji — h + hh ipi — ai + a^di^F 

J2 ^ h + hi?, ■02 = ^2 + a^di^F 

J3 = F{h,l2) + hh = l-a - Lut, 

and H = H' + lu I3. Thus, near the resonant phase-space surface, J = Jr + 
and we can expand the new Hamiltonian around this surface, 

H'ixPJ) = HoiJR + I3I) + V{Jr + hi) cos^3 - huj (2.67) 

= HoiJR.) + h il-n{jR) -u) + ^Dli 

+ ViJR + hi) cos4>3, (2.68) 

with D{Jr) = khdj,dj^Ho{JR). 

If the variations of H{J) are small around Jr, we can replace V{Jr + h n) by 
V{Jr) in the latter expression, and obtain the new Hamiltonian H' = Hq{Jr) + 
Hi..jR{i^3,h), with 

Hij^i^,!) = ]^DI^ + Vcosxl,. (2.69) 

Thus, the problem has been reduced to a ID Hamiltonian problem for the angle-action 
variables (-0, I)={'4>3, h)- 



Substituting the expression of the TAE perturbation, Eq. (2.62) into Eq. (2.571, 
we obtain 

iP = pO + nC ~ Lut, (2.70) 

J ^ JCZICR ^ (2.71) 
?i n 

D « n'^iJn), (2.72) 

V = -ieCit)Vp{JeR,Jcn), (2.73) 
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where the subscript R means that the quantity is evaluated at the resonance, and 
with a = 9, C), J = (J^, Jo, Ji^), and I — (0,p, n). 



2.3.2 Reduction of the collision, source and sink terms 

A first, simple model is obtained by reducing the effects of collisions to the recovery 
of an equilibrium energetic particle distribution, with a recovery rate Vaiv)- 

A more rigorous treatment of collision processes is obtained if we project a collision 
operator that describes Coulomb collisions perceived by energetic ions, on the resonant 
phase-space surface. We consider collisions on energetic particles by thermal electrons 
(s = e), ions (s = i), and carbon impurities (s = c), and describe them by a Fokker- 
Planck collision operator [HS112j that acts on the distribution of energetic particles 
(s ~ b). In spherical coordinates (f,0), where 9 is the angle between v and b, 
neglecting gyroangle dependency, 



dt 



1 1 



t'dcfl 



d 



coll. 



2 sin 6 de 



sin O 



dl 



dv 



1 df 



(2.74) 

where i/defl, ^^siow and are pitch-angle scattering, slowing-down, and parallel velocity 
diffusion rates, respectively, = v cos Q is the parallel velocity of energetic particles. 

We consider a TAE with toroidal mode number n, resulting from the coupling of m 
and m + l poloidal modes. To simplify the following discussion, we consider strongly 
co-passing beam particles that resonate with the TAE at a velocity v « v\\ = va- 
Then the resonance condition is given by 57 = wa, where fl is given by Eq. (2.6). To 



project the Fokker-Planck operator on the resonant phase-space surface, we follow the 
procedure described in Refs. |BBP97a1 ILBS09| . We replace d^f by Jbodnf, where 
J is the Jacobian of the coordinate transformation from v\\ to 51, 



J 



2r2ebBo 



(2.75) 



and S is the magnetic shear. Here Cs and nis are charge and mass of a species s, 
respectively, and b stands for beam particles. This procedure yields 



dt 



coll. 



2 df 



(2.76) 



with 



^slow 



j2 {v« cose 



^'dcfl) 



z^dcfi sm 



e) 



(2.77) 
(2.78) 



We assume Maxwellian background distributions, with a same temperature Tq. 
Typical experiments satisfy the following ordering of thermal velocities, vtc < VTi ^ 

^ VTe, while the beam energy Ef, is much larger than Tq. Within these assump- 
tions, around the resonance. 



VllJ 



2v^ 



E 

s 

E 



erf 



2fh 
2Vs 



2vu] erf 



(2.79) 



(2.80) 
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where rjs = v/vts, = vsinQ, v\\ = va, 



lbs 



elel log A 



(2.81) 



eo is the vacuum permittivity, and log A is the Coulomb logarithm. Since the magnetic 
moment is an invariant of the motion of injected beam ions from deposition to resonant 
phase-space surface, = u^(l — R^/FLl), where Vb is the velocity of beam particles, 
and Rt is the tangential radius of the beam. 

The equivalent collision operator in the BB model is obtained by substituting 
r2 = few in Eq. (pjel). 



2.3.3 Reduction of the background damping mechanisms 

Since the time-scale of fast-particle evolution is much faster than background thermal 
populations evolution, these two dynamics are decoupled. Hence we can reasonably 
treat the effects of background damping in an extrinsic way. We assume that all 
background damping mechanisms affect linearly the wave energy W, 

dW 

— = -2-idW{t). (2.82) 

Damping includes mechanisms such as radiative damping, which strength depends 
on the frequency [MM92 . Hence, in a rigorous model, 7^ should be a function of w. 
However, theory needs to be developed before this complex interplay can be taken 
into account. Thus, we limit our framework to cases where 7^ can be treated as a 
constant. This framework is consistent with a fixed mode frequency. 



2.3.4 Limitations 

We assumed an isolated single mode, which is reasonable for sufficiently core-localized, 
low-n TAEs. However, for future devices with higher-n, an other model that includes 
multiple-modes interactions has to be developed. 

Since it assumes a fixed mode structure, implying a fixed Magneto-Hydrodynamic 
(MHD) equilibrium, the above reduced model looses its validity on a time-scale of 
MHD equilibrium evolution, which is of the order of the second on large devices 
such as JT-60U. We must also require that wave amplitude is low enough, such that 
nonlinear redistribution of energetic particles does not significantly alter the mode 
structure and frequency. In practice, since the frequency and growth rate of TAEs 
are very sensitive to the q profile, if the frequency of a low-amplitude TAE observed 
in experiments stays nearly constant during a certain time-window, we infer that the 
fixed-niode-structure assumption is satisfied for this time- window. 

In the case of frequency sweeping, which is the case we apply to experiments, it is 
sometimes argued that since the frequency is changed, so must be the mode structure. 
However, we must distinguish at least three classes of frequency sweeping, namely, 
slow frequency sweeping (slow-FS), fast frequency sweeping (fast-FS), and so-called 
abrupt large-amplitude events (ALE) [SKT+Ol] , although it is not clear for the latter 
if the frequency does sweep. In the case of JT-60U shot number E32359, which is 
analyzed in Chap. [5] slow-FS have a timescale of 100 — 200 ms, and their frequency 
is correlated with bulk equilibrium variations, therefore they are out of the scope of 
the above reduce model. Fast-FS have a timescale of 1 — 5 ms, and the associated 
redistribution of energetic ions is relatively small |STI+02) . therefore are consistent 
with a fixed-mode-structure hypothesis. Although the occurrence of fast-FS and 
ALEs seems to be linked, ALEs are identified as so-called Energetic Particle-driven 
Modes (EPMs) [BFV+Of] . have significantly larger amplitude and shorter timescale 
(200 — 400 /is), and induce significant loss of energetic ions, and are out of the scope 
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of this work since we assume a weak drive and a constant density of energetic ions in 
the BB model. 

Finally, modeling all background damping mechanisms as an extrinsic, fixed linear 

damping on the wave is a strong assumption, whose validation requires more under- 
standing of these mechanisms. We must assume that docs not depend neither on 
the wave amplitude, nor on the energetic population. In the case of frequency sweep- 
ing, the assumption is clearly violated if the nonlinear modification of frequency is 
of the order of the linear frequency, especially if a chirping phase-space structure ap- 
proaches the SAW continuum, where damping rate depends largely on the frequency. 

Overall, the above reduced model is suitable for describing resonant interactions 
between energetic particles and a weakly driven, isolated TAE, even for slightly- 
chirping modes, as long as 

• phase-space structures are well confined within the continuum gap ; 

• redistribution of energetic population is negligible as far as wave dispersiveness 
and damping mechanisms are concerned ; 

• we look at time-scales much smaller than the equilibrium evolution time-scale. 
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Chapter 3 

The Berk-Breizman model 



The BB model, describing the fully nonlinear wave-particle interactions between ener- 
getic particles and an electrostatic wave in a ID plasma, is an extension of the Vlasov- 
Poisson system, where we take into account collisions and external wave damping. In 
this system, we consider a bump-on-tail velocity distribution comprising a Maxwellian 
bulk and a beam of energetic particles, and we apply a small electrostatic perturba- 
tion. The apparent simplicity of the corresponding equation system hides surprisingly 
rich physics. Depending on the parameters of the model, the perturbation may be 
damped or amplified due to a transfer of energy between resonant particles and the 
wave. In the stable case, when the perturbation is small, linear theory predicts expo- 
nential decay of the wave amplitude, which in the absence of collisions and external 
damping is known as Landau damping jLan46| . In the unstable case, on the contrary, 
linear theory predicts exponential growth of the wave amplitude. Then, trapping 
of resonant particles significantly modifies the distribution function and an island 
structure appears. Saturation of the instability and following nonlinear evolution 
are determined by a competition among the drive by resonant particles, the external 
damping, the particle relaxation which tends to recover the initial positive slope in 
the distribution function, and particle trapping that tends to smooth it. It has been 
predicted jBBP97a[ lBBP+97c[ IBBP96j and observed |VDR+03| that three kinds of 
behavior emerge, namely steady-state, periodic, or chaotic responses, depending on 
the strength of each factor. In addition, chaotic solutions can display significant shift- 
ing of the mode frequency (chirping) , both upwardly and downwardly, as holes and 
clumps are formed in the distribution |BBP97b| IBBP98| lBBC+99] . 

In Sec. |3.1[ we recall the equations of the BB model, in both full-/ and 6f ap- 



proaches. In Sec. 3.2 we show the linear dispersion relation, and present tools that we 
use for accurate linear analysis. For our purposes of validating and extending BB the- 
ory, and of applying it to TAEs, we develop a kinetic code based on the Constrained 
Interpolation Profile - Conservative Semi-Lagrangian (CIP-CSL) scheme INTYTOl] 
for solving the initial value problem. In Sec. |3.3( we present the main principles of our 
code, which we name COBBLES. In the full-/ case, we show that a locally conservative 
implementation is a key point for robust simulations in experimentally-relevant con- 



ditions, which are particularly stringent in a numerical point-of-view. In Sec. 3.4 we 
verify nonlinear capabilities of COBBLES. In the coUisionless limit without external 
damping, we are able to solve the simpler Vlasov-Poisson system and recover satura- 
tion level, relative oscillation amplitude, and the so-called Bernstein-Greene-Kruskal 
(BGK) steady-state solution (BGK5 7> Then, we analyze the conservation and con- 
vergence properties for a system with finite 7^ and collision frequencies. Further, 
we benchmark COBBLES against a parameter scan given in former works by Vann 
[VDR"'"d3] . A drag and diffusion collision operator is verified by recovering qualitative 
steady-state distributions predicted by theory. Finally, we consider multiple- waves in- 
teraction and estimate a particle diffusion coefficient which agrees with quasi-linear 
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Tab. 3.1: Correspondence between physical and normalized quantities. 



Physical quantity Normalization constant 

Time lo^ ^ 

Length Aq 

Velocity Vth 

Density tlq 

Distribution / no/vth 

Electric field mUj^/(eAD) 

Energy, Hamiltonian "^^t/i 

Power mX-aTiQvjf^uip 



theory. In Sec. 3.5 we summarize the analogies between BB model and ID model of 
TAE. 



3.1 Basic equations and physics 

Depending on the application, it may be convenient to cast the BB model either in a 
self-consistent form (full-/) or in a perturbative form (<5/). 

3.1.1 Normalization 

For the sake of concision in this thesis, and to avoid numerical treatment of too large 



and too small numbers in simulations, we adopt the normalization shown in Tab. 3.1 
where the plasma frequency is defined by uf^ = rtoe^/(eoTO), e, m, and no are the 
charge, mass, and total density, respectively, of the evolving species. Ad — Vth/^p is 
the Debye length, and is a typical thermal velocity. 



3.1.2 Full-/ BB model 

We consider a ID plasma with a distribution function f{x,v,t). In the initial condi- 
tion, the velocity distribution, 

foiv) s 7(«,t = 0) - f^'iv) + fo^iv), (3.1) 
where / is the spatial average of /, comprises a Maxwellian bulk, 

fo'Hv) - ^^e-K^)\ (3.2) 
vtm V 27r 

and a beam of high-energy particles, 

foi^) = ^^e-K^)', (3.3) 
vtb V ^tt 

where um and ub are bulk and beam densities, which verify um + ns = 1, vtm 
and Vtb are thermal velocities of bulk and beam particles, and is the beam drift 
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velocity. To ensure charge neutrality, we assume a fixed background population of 



the opposite charge with a distribution function fo{v). Fig. 3.1 shows two typical 
initial distribution functions, with a cold bulk and a weak, warm beam. The first 
one, to which we refer in the following as distribution A, with tib — 0.1, vtm = 0.5, 



vtb = 1-0, and vb = 4.5, is used to benchmark our code in 3.4.3 The second one, 
distribution B, with ub — 0.1, vtm = 0.2, vtb = 3.0, and vb — 5.0, hence a warmer 
beam and a colder bulk, is used to validate and develop some aspects of BB theory. 
For both distributions, we will always choose a system size L ~ 2t: /k with k = 0.3. 

The evolution of the distribution is given by the kinetic equation 

at +"ax+^a^ =c(/-/o), (3.4) 

where E is the electric field, and C{f — fo) is a collision operator. 

In this work, we consider either one of the following two collision models. On the 
one hand, a large part of existing theory for the BB-model deals with collisions in the 
form of a Krook operator jBGK54j , 

Ck(/-/o) - -'^aiv) if - h) , (3.5) 

which is a simple model for collisional processes that tend to recover the initial dis- 
tribution at a rate i^^, including both source and sink of energetic particles. If we 
assume cold and adiabatic bulk plasmas, Vaiv) acts only on the beam. Reflecting this 
situation, we design the velocity dependency of Ua (v) as 

else ' (^-^^ 

where satisfies fo\viy)/fo''{0) = e^, and we choose — 10^^. Hence, i^a{v) is 
constant in the beam region, and zero in the bulk region, except for the benchmark 



in Sec. 3.4.3 where it is explicitly stated to be constant everywhere. 
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On the other hand, a more reaUstic colhsion operator, the one-dimensional pro- 



jection of a Fokker-Plank operator, Eq. (2.761, includes a dynamical friction (drag) 
term and a velocity-space diffusion term. 



Cfp(/-/o) 



'^/H5(/-/o) 



k 



dv 



fc2 dv'^ ' 



(3.7) 



where k is the wave number for the resonance under investigation, and with similar 
velocity-dependence for Vf and j/^j. An other large part of existing theory deals with 
the latter operator in the absence of drag {vf = 0). Investigations of the effects of 
dynamical friction are fairly recent [LBSOQl ILFCIO) . 

We define the effective collision frequency as = Va in the Krook case and 
i/cS = i^d/lLo i'^ the case with diffusion. A dimensional analysis gives a typical 
lifetime of phase-space structures as . 

In the expression of the electric field, 



E{x,t) = Ekit)e 



tkx 



+ C.C.. 



(3.8) 



we assume a single mode of wave number fc, reflecting the situation of an isolated 
single mode AE. The displacement current equation (DCE), 



dE 



= - J v{f-fo)dv - 2jdE, 



(3.9) 



yields the time evolution of the wave. In the initial condition we apply a small per- 
turbation, f(x,v,t = 0) = /o(w)(l -l- ecosfcx), and the initial electric field is given by 
solving Poisson's equation. In Eq. (3.9), an external wave damping has been added to 



model all linear dissipation mechanisms of the wave energy to the background plasma 
that are not included in the previous equations [BBY93j. The presence of a factor 2 
in front of 7^ is consistent with Berk and Breizman's literature and will be justified 
in Sec. [X2I 



Conservation properties 



Before deriving the conservation properties of this model, it is useful to note the 
following property. If f{x,v,t), g[f{x,v,t),t] and h{x,v,t) are arbitrary functions, 
analytic in a phase-space F = [x, v), then 



da; 



djgh) 
dx 



dv 



djgh) 

dv 



= 



for usual boundary conditions. 



(3.10) 
(3.11) 



where integration in the l.h.s. is over the whole phase-space surface. 

In the ideal situation, the model presented above ensures conservation of total 
particle number N(t) = J fdxdv. This is proven by taking the integral over the 
whole phase-space of the kinetic equation, which can be written as 



dtf ~ = C(/-/o), 

where h is the Hamiltonian, 

h{x,v,t) = u^/2 + ip{x,t), 
and if is the electrostatic potential. Thus we obtain 



diV 
~dt 



J C{f-fo)dxdv. 



(3.12) 



(3.13) 



(3.14) 
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With Krook operator, if Vai^) = Va is taken as constant, the latter equation can be 
written as 

+ Va^N = 0, (3.15) 



where AA^(i) = N{t) ~ iV(0). Since AA^(O) 0, Eq. ( |3.15[ ) yields the conservation 
of total particle number, dAA^/di = 0. However, in numerical simulations, some 
spurious leakage of particles from velocity boundaries induces a small error in this 
conservation. When v^iv) has the velocity dependence of Eq. (3.6), Eq. (3.15) is 
changed to 



dt 



(/ - /o) Av. 



(3.16) 



In the bulk part u < Wj/, we assume that the variation of the distribution is negligible, 
|/(w,t) — /o('^)| ^ fo{v), and we show the approximative conservation of total particle 
number. 



dAA^ 



dt 



(3.17) 



With Fokker-Plank operator, dN/dt = is immediate from Eq. (3.14). 

Let us now derive an energy equation to relate power transfers between wave, 
particles, and external damping. Detailing each term is useful to clarify different 
possible decomposition of the power transfer, corresponding to different point of view. 
Taking the integral over phase space of the product of the Hamiltonian with the kinetic 
equation, and dividing it by the system size yields 



1 
L 



h^dV 
dt 



v"^ df dx 



dv 



df dx , 



P. 



2 dt L 

where the right hand side shows the coUisional power transfer = Pj + P^ , 



C(/-/o)y^d., 



C(/-/o)<y5^du. 



(3.18) 

(3.19) 
(3.20) 



The left integral of the left hand side is the kinetic part dT/di, where T is the total 
particle kinetic energy density, 



m 



2 



(3.21) 



Substituting the kinetic equation into the right integral of the left hand side, we find 
out that the field part is —Ph + P^f, where we define the particle power transfer as 



Ph{t) 



V E f — dv. 



(3.22) 



There are two ways of decomposing Ph. The substitution of DCE (3.9) into the 
expression of Ph yields 

Ph = -^-Pd. (3-23) 
where the electric field energy density is defined by 

and the external damping power transfer by 



dx 



(3.24) 



Pd{t) = 47d£. 



(3.25) 
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On the other hand, by substituting the kinetic equation (3.4 1 into the expression of 
dT /At we obtain 

Ph - -Pj + (3.26) 

We can further separate the response of the particules into resonant and non- 
resonant pieces if we decompose the distribution into / = + /^^ and define 
T ^ + T^^. The non-resonant part of the kinetic energy is the sloshing energy, 

-^NR ^ j Y^-jm^^ (3.27) 

For non-resonant particles, the velocity is oscillatory and we can replace the amplitude 
of its oscillation by the linear response Eq/uj, where Eq = 2|i?fc|, and we obtain 
-y-NR _ g rpj^g wave energy W is composed of the field energy £ and the sloshing 
energy 7~^^ which supports the wave, 

W = + £ = 2£. (3.28) 

Finally, we can express the power transfer equation in two equivalent ways, de- 
pending on whether we consider the non-resonant kinetic energy as part of the total 
kinetic energy or as part of the total wave energy. In the electric field point of view, 

— + + 47rff = 0, (3.29) 
at 

showing the balance between the field and all the particles. In the wave-as-quasi- 
particles point of view, 

^ + P« + 27rfW = 0, (3.30) 
where Pj^ is the resonant power transfer, 

^ JvEf'fdv = P, (3.31) 

showing the balance between wave and resonant particles. 

3.1.3 6f BB model 

If the bulk particles interact adiabatically with the wave, their contribution to the 
Lagrangian can be expressed as part of the electric field. Then it is possible to adopt 
a perturbative approach, and to cast the BB model in a reduced form that describes 
the time evolution of beam particles only jBBP95a[ ICD93| . The evolution of the beam 
distribution, f^[x,v,t), is given by the kinetic equation 

where the pseudo-electric field E is defined as 

E{x,t) = Q{t) cos{ip) - P{t) sin(V'), (3.33) 

where ip = kx — ut. In this model, the real frequency of the wave is imposed as 
w = 1. This restriction does not forbid nonlinear phenomena like frequency sweeping, 
since both amplitude and phase of the wave are time-dependent. In this thesis, we 
renormalize physical quantities for the Sf model so that they do not depend on k. 
In practice, we choose fc = 1. In the collision operators, i^f and vj^ are taken as 
constants, since, with the 5f description, velocity dependency is not needed to avoid 
affecting bulk plasma with collisions. 
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(V-Vr) / Vp 



Fig. 3.2: Initial velocity distribution function, normalized to /jv = 7lo/(v Svrufl). 



The evolution of the pseudo-electric field is given by 

J f^{x,v,t) cos{^)dxdv ~jdQ, (3.34) 

^ j f^{x,v,t) sin(V;) dx dv - 7^ P. (3.35) 

The initial values of Q and P are given by solving Poisson's equation. Note that the 



dQ 
dt 
dP 



latter equations, without factor 2 in front of 7^, are consistent with Eq. (3.9) 



In the coUisionless case, one can see from the linear dispersion relation Eq. (3.58) 
that = 1 only if /(f is symmetric around the resonant velocity, vr = Lo/k. Since we 
assumed w = 1 from the start, we consider only such distributions, for the model to 
be self-consistent. The velocity distribution of beam particles in the initial condition 
is shown in Fig. |3.2[ A constant slope is imposed between v = —v^ and v = v^,, where 
Vc is an arbitrary cut-off velocity. The zero average ensures that the plasma frequency 
is not perturbed by the beam density. Smooth joins between the constant gradient 
region and the large velocity regions are necessary to prevent numerical oscillations at 
V » ivc- Since we always choose Vc large enough so that border effects are negligible, 
an initial distribution is fully characterized by its slope at resonant velocity, in other 
words by 7z,o- 



Conservation properties 

Arguments similar to those for full-/ model yield the conservation of total particle 
number. The power balance is changed to, 

Ph + Pe + Pd ^ 0, (3.36) 
where Ph is the kinetic power transfer, 

Ph ^ j vEf^ dxdv, (3.37) 

Pe is the electric field power transfer. 
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Tab. 3.2: Non exhaustive list of approaches and coUision operators in the hterature. 
"BB" refers to Berk, Breizman and coworkers, "Lesur" refers to this thesis. 



Author Approach 



Cohisions 



BB Sf Krook / Diffusion 

Vann full-/ Krook 

Lilley Sf Krook / Diffusion / Diffusion+Drag 

Lesur Sf / full-/ Krook / Diffusion / Diffusion+Drag 



and Pd is the power transfer due to external damping and collisions. In the Krook 
case, 



Pd = 



27r 

T 



(7d 



PQ - QP) + 7d (0' + P^ 



(3.39) 



while in the Fokker-Planck case, collisions do not contribute to this latter power 
transfer, thus Pd is obtained by substituting j/q = in Eq. (3.39). 

Compared to the full-/ model, the Sf model does not take into account effects of 
time-evolution of bulk particles, which is a caveat when assessing limit of theory that 
breaks-up when phase-space structures approach the bulk, but it has an advantage 
in the application to experiment, where we assume fixed mode structure, hence fixed 
background plasma. Moreover, the velocity range required to simulate a similar reso- 
nant region can be significantly reduced with the Sf model, saving computation time. 
Since, in this thesis, we often refer to literature by Berk and Breizman, by Vann, or 
by Lilley, we clarify which approaches and which collision operators have been studied 
by these authors, in Tab. 3.2 



3.2 Linear analysis 

When the perturbation is small, linear theory predicts exponential growth or decay 
|Lan46| of the wave amplitude. For the full-/ model with Krook collisions, the linear 
dispersion relation, 

7 + 27. - .a; = / , ^ ^f^'^" .dv, (3.40) 

where T is the appropriate Landau contour |Lan46| . yields the linear growth rate 7, 
and the real frequency uj of the wave. We implemented an algorithm to solve the 
latter equation, applying a method of residue for locating the zeros of an analytic 
function in the complex plane |Dav86| . We refer to this algorithm as Davies solver. 
In the following, jl is defined as the linear growth rate for 7. = = 0. 

In the coUisionless limit, if we assume a small perturbation and a linear growth 
rate 7 much smaller than the real frequency u, the dispersion relation reduces to 

7 = 7-L ^| , (3-41) 



where 



and V is a notation for Cauchy's principal value. In the cold Maxwellian limit, 
7i — 7lo, where 7^0 is a measure of the slope of initial distribution at resonant 
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Fig. 3.3: Growth rate in the cohisionless hmit. Sohd hne corresponds to Eq. (3.44 1, 
dashed hne corresponds to distribution A, for which 7^ = 0.1981, and dotted hne 
corresponds to distribution B, for which 7^ = 0.0324. 



vefocity, 

_ TT dfo 
" 2fc2 dv 



(3.43) 



In this hmit, a simple relation, 7 = 7l — 7d, stands. However, with our choice of 
distribution function, we must keep in mind that there is some discrepancy between 
7 and 7l - Jd, 

7 = 7L f 1 - —] for i^a = 0,l « ^- (3.44) 
Fig. |3.3| shows the growth rate estimated by Davies solver as a function of external 



damping, in the cohisionless limit, for both initial distributions A and B. Eq. (3.44) 
is recovered in the limit of small 7. 

With Fokker-Planck collisions, the kinetic equation in Fourier-Laplace space is a 
second order differential equation in v, which prevents a similar treatment. However, 
we can take another approach, which is also valid in the Krook case, where we search 
for solutions of the form exp(pt), where p = -f ~ luj. Writing fk{v, t) — fp{v)d'^ and 
Ek{t) = EpcP* the Fourier component of f — fo and E, respectively, we obtain a linear 
equation system, 

^P + 'k-)fp + Ep— ^ -.,/^ + __ + __, (3.45) 
{p + 2-ia)Ep = -Jvfpdv. (3.46) 

Discretizing the velocity space as vj = jAv for j = 1---Ni,, the latter system is 
approximated to first order accuracy in Aw by 

dfo 

{p + lkVj)fj + Ep -T^iVj) = -Vafj + {fj + l - fj-l) 

+ ^2^^ -2/.+ /.-i), (3.47) 

iV„ 

{p + 2jd)Ep = -Avj^vjfj: (3-48) 



28 



Krook 






(a) 






X 








X 






1 


128 X 




-2 -1 





1 


2 




0) 







0.2 



-0.2 
-0.4 
-0.6 
-0.8 
-1 
-1.2 
-1.4 



Fokker-Planck 

5K 



(b) 




/ 



ypy!j ! ii j. i i(i ii p iiiii 4. i '" ' mHHttitiiffi 



Nv=128 
Nv=256 



-2 



-1 





CO 



Fig. 3.4: Eigenvalues for initial distribution A, with 7^ = 0.1. (a) Krook case, with 
I'a = 0.25. (b) Fokker-Planck case, with i/f = 0.02 and I'd = 0.05. 



where fj = fp{vj), and boundary conditions are /o = Jn^.+i = 0. This system of 
Ny + 1 equations can be put in matrix form, 

M • F = pF, (3.49) 

where F is a Ny + 1 dimension vector defined by 

F, = f, {j^l---Ny), (3.50) 
Fjv„+i = Ep, (3.51) 

and M is a Ny + 1 dimension square matrix defined by 

M,j = -^kvj - ua - {j^l...Ny), (3.52) 

= ^ + ^ U-l-Ny^l), (3.53) 

M»+i - pI;^ - ^ = (3.54) 
MAr,.+ij = -9./o(^'j) (i = l---^„), (3.55) 
Mj,JV„+i = -I'jAv (j = 1...7V,), (3.56) 

MAr„+i,Ar„+i = -2jd, (3.57) 

where M, j is the element of column i, line j. We solve the above eigenvalue problem 
using LAPACK library. In the Krook case, w + 17 = few — iVa constitutes a continuum 
of trivial solutions with E — 0. This eigenvalue method does not yield solutions with 
7 < —I'a- Moreover, as 7 approaches —i'a, increasing number of grid points are needed 
to accurately estimate the growth rate, since continuum solutions tend to perturb 



nontrivial solutions. Eigenvalues found by this method are shown in Fig. 3.4 'a), for 
initial distribution A, for which 7^ = 0.1981, and with — 0.1 and i^a = 0.25. The 
least stable solution is w = 0.9953, 7 ~ —0.04551, which is confirmed with Davies 



solver (see Table 3.3). Thus, damped solutions can be found, if I7I < Va- In the 
Fokker-Planck case, the continuum is topologically changed, and depends on both 
distribution function and grid points number. Thus, when estimating growth rates 
with this method, we must be careful to use values that are converged with Ny. This 
is illustrated in Fig.|3.4||b), which shows eigenvalues for distribution A, with "fd — 0.1, 
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Vf = 0.02 and = 0.05. The most unstable solution is lo 
which is in agreement with full-/ simulation (see Table 3.3 1. 



1.0204, 7 = 0.1737, 



For the 5f model with Krook collisions, the linear dispersion relation is changed 



to 



7 + 7d - « (t^ - 1) = 



1 

2k 



dv. 



(3.58) 



/r (7 + i'a) + « [kv - uj) 

As announced before, we have lu — 1 only if the imaginary part of the right-hand side 
vanishes, in other words if is anti-symmetric around the resonant velocity. In the 
collisionless case, for 7/ti; <C 1, Eq. (3.58) yields 



7 



ILO 



Id- 



(3.59) 



If we search for solutions of the form exp(pt), as fk{v,t) = /p(u)eP* and Z{t) 
Zpe^*, where Z{t) = [Q{t) + tP{t)] exp(— it), we obtain a linear equation system. 



Zpdfo 
2 dv 



{p + ikv)fp 

{p + 27d 4- t)Ep 



-Va fp 
1 



_dfp 
k dv 



fc2 dv"^ 



fpdv. 



(3.60) 
(3.61) 



The discretized version of the latter system can easily be put in the form of an 
eigenvalue matrix problem, and solved in a way similar to the full-/ case. 



3.3 The kinetic code COBBLES 
3.3.1 Numerical implementation 

Let us recall that the BB model is an extension of the Vlasov-Poisson system, which 
is recovered in the collisionless, closed system (7^ = 0) limit. In a previous work 
|LIT07| . we developed a ID semi-Lagrangian full-/ Vlasov code, based on the Cubic- 
Interpolated-Propagation (CIP) scheme |NY99| and the splitting method |CK76| . 
which enabled accurate simulations of the Vlasov-Poisson system. In this thesis, 
we extend our code to include distribution relaxation and extrinsic dissipation, and 
develop a 6f version. We refer to these codes as full-/ COBBLES and 6f COBBLES, 
respectively, COBBLES standing for Conservative Berk-Brcizman semi-Lagrangian 
Extended Solver. 

In both codes we solve DCE instead of Poisson equation. Looking at the spatial 



average of Eq. (3.9 1, 



dE 

lit 



= - J v{f - fo)dv - 2jdE, (3.62) 



we see that a small deviation from a constant total momentum can be the source 
of a systematic error in the average electric field. Such deviation arises when Krook 
collisions are included, or can be caused by numerical error. To avoid this problem. 



we replace / vfodv by / vfdv in the DCE |Van02] . Then Eq. (3.62) is changed to 
dtE + 2jdE — 0, which ensures a zero average electric field, since E\^_^ — 0. 

Let us now describe the main points of our algorithm. All quantities like / are sam- 
pled on uniform Eulerian grids with and Ny grid points in the x and v directions, 
respectively, within the computational domain {{x, v) \ < x < L, v^^ £ w < Wmax}- 
For distribution A, cut-off velocities are always chosen as Wmin — —8 and Wmax — 8. 
For distribution B, Wmin = —10 and Wmax — 18. We define the Courant-Friedrichs- 
Lewy number CFL = fmax At iV^, / (2L) as a measure of the time-step width At. We 
use the Strang splitting |Str68j formula to obtain a second-order accuracy in time 
jVDR"'"03 . For each time-step, we perform the following steps, 
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Fig. 3.5: Strong-scaling for the COBBLES code, with and without OpenMP enabled. 



1. Advect dtf 
2 



vdxf = for a time At/ 2 



Solve dtf = -ya{f-h) - {y}/k)d,fo + {vl/e)dl{f - fo) iov ixiiuve 
M/2 



for a time At 



Solve DCE for a time At/2 
Advect dtf + {E{x) - v}lk)d^f 
Repeat the step 3. 

6. Repeat the step 2. 

7. Repeat the step 1. 

Numerically, step |3] is performed by a forward Euler scheme. Note that the 
implementation of friction is split into steps [2] and |4j In step [2j / is replaced by 
/o + exp(— i^aAt/2)(/ — /o) — {i'jJk)dyfoAt/2, then the diffusion equation is solved 
by the Crank-Nicolson method |CN47j . The remaining problem, corresponding to 
steps 1 and 4, is the advection of a ID hyperbolic equation, 



dtF + udxF 



0, 



(3.63) 



where u is constant in the A direction, A is a generalized advection coordinate, and 
is a general function of A and t. We aim at long-time accurate simulations in the 
whole (7^, Va) space. The choice of advection scheme is of crucial importance to 
reach this goal. In Appendix [C] we recall the CIP-CSL algorithm, which we use for 
solving Eq.(|3.63 l, and its extension to the position- velocity phase-space, as presented 
in Ref. |NTYTOl] . The key idea is that in addition to the distribution function, 
we advect its integrated value p to keep a flux balance between neighboring grids. 
We justify this choice in the following section. Boundary conditions are periodic in 
configuration space, and zero-flux at velocity boundaries. 

COBBLES is coded in Fortran 90 language. It is parallelized in a hybrid fash- 



ion, using MPI in the velocity direction and OpenMP. Fig. 3.5 shows the speed-up 
on JAEA's BX900 systems as the number of Central Processing Units (CPUs) is in- 
creased, at fixed grid-points number x Ny = 6A x 4096. Although there is room 
for optimization, the observed scaling properties are sufficient for our purposes. The 
use of OpenMP in addition to MPI provides a significant speed-up for 128 CPUs. 
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Differences between 6f and full-/ versions are in the init ializa tion, which defines 
instead of /, and in the DCE, which is replaced by Eqs. (3.34) and (3.35). 
Note that we also implemented the treatment of two species, which was used 

to benchmark simulations of Ref. [NLG+IO] . though this feature is not used in this 

thesis. 



3.3.2 Comparison of advection schemes 

Compared to the benchmark of |3.4.3[ our choices of parameters for testing theory, 
in Chap. |4] constitute difficult conditions for numerical stability of full-/ simulations 
and drastically increases computational cost. As a consequence, we must take special 
care in choosing an advection scheme that minimizes computational time. Therefore 
we discuss the relevancy of the CIP-CSL advection scheme, and present a comparison 
with four other advection schemes. 

In choosing the advection scheme, we focus on stability and convergence proper- 
ties, which are estimated with severe benchmark parameters relevant for analysis in 
Sec. 4.1.2[ where we use distribution B (with a cold bulk and a weak, warm beam). 



Compared to distribution A, which is used as initial condition in the following bench- 



mark (Sec. 3.4.3), simulations with initial distribution B are more sensitive to numer- 
ical errors such as numerical diffusion. Firstly, the colder the bulk, the less grid points 
are available in the bulk, leading to artificial heating. Secondly, the weak warm beam 
induces weaker linear instabilities, which produce narrower islands in phase space. To 
resolve such a narrow island, increased grid resolution is needed. Furthermore, for 
steady-state solutions, when the island is narrower we observe unphysical drive after 
nonlinear saturation, which suggests that the region of flattening acquire spurious 
gradient by influence of surrounding distribution. In this work, we aim at producing 
a numerical scan of nonlinear behavior in the whole parameter space. Near marginal 
stability, the linear growth rate 7 is very small (we limit the investigation range to 
I7I > 10~^ to avoid excessive computation cost) and long-time computations {t ^ 10^) 
are required. For this reason, we cannot afford too much grid points, and we have to 
take utmost care in choosing a robust and quickly converging numerical scheme. 



A comparison of several advection schemes for one of the case of Fig. 4.2 (with 



distribution B) is shown in Fig. 3.6 The time evolution of a beam instability with 
a low dissipation i/aiv > v,y) = 0.002, and a small external damping — 0.002, for 
increasing grid resolution, is compared to a reference run for each of five schemes : 
Flux-Balance (FB) [Fij99l, CIP |NY99] . CIP with rational function interpolation (R- 
CIP) LXYNI96 , CIP-CSL, and Rational - CIP-CSL (R-CIP-CSL) scheme jXYPK02] . 
The reference run is obtained with a high resolution x — 256 x 4096 using 
CIP-CSL. The CIP scheme is a low-diffusion and stable scheme, and is implemented 
in a way that exactly conserves the total mass. However, it is not locally conservative. 
After several amplitude oscillations in the nonlinear phase, we observe the apparition 
of numerical oscillations in the velocity direction in a large gradient region of the 
distribution, which appears between a cold bulk and a beam. While, in this test case, 
numerical divergence eventually occurs even for very high resolution with the CIP 
scheme, the other schemes show convergence to a same solution. The FB scheme 
is only second-order accurate, so that convergence is slow compared to the CIP- 
based schemes, which are third-order accurate in general |NY99| . Rational function 
interpolation aims at preventing numerical oscillations by preserving convex-concave 
and monotonic properties, but at the expense of this property, numerical diffusion 
produces spurious drive leading to higher saturation levels. R-CIP-CSL produces 
less numerical diffusion than R-CIP, but convergence is slower than with CIP-CSL. 
Finally, the CIP-CSL scheme shows quick convergence without unfavorable numerical 
oscillations, and therefore, we use this scheme in the following simulations. 

As an illustration of 5f model, we include in Fig. |3.6|the time-evolution of beam 
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Fig. 3.6: Time evolution of the normalized bounce frequency for different advection 
schemes. Solutions are shown for grid resolution x of 64 x 512, 128 x 1024, and 
128 X 2048, and for a reference run described in the text, (a-e) Full-/ simulations with 
initial distribution B, Va{y > Vu) = 0.002, 7^ = 0.002, and CFL ^ 0.9. Each sub- 
figure corresponds to one of five advection schemes, (f ) 5f simulation with 7^0 = 0.037 
(such that 7L is the same as for full-/), = 0.002, 7^ = 0.002, and CFL = 0.9. 
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Tab. 3.3: Linear frequency and growth rate, obtained by solving the eigenvalue prob- 
lem, or using Davies solver, or by fitting an exponential and performing a Fourier 
analysis of electric field time-series in full-/ COBBLES simulations. In both cases, 
the initial distribution is A, with 7^ = 0.1981, and 7^ = 0.1. 



Collision operator 


Krook 




Diffusion-KDrag 


Collision frequencies 


Ua = 0.25 




Vf = 0.02, Vd = 0.05 


Eigenvalue 


a; = 0.9953, 7 = - 


-0.04551 


UJ = 1.0204, 7 0.1737 


Davies solver 


UJ = 0.9953, 7 = - 


-0.04551 




COBBLES simulation 


w = 0.9948, 7 = - 


-0.04549 


UJ = 1.0176, 7 = 0.1743 



instability obtained by 5f COBBLES (with the CIP-CSL scheme), for similar param- 
eters as those of the full-/ simulations. 

3.4 Verification of COBBLES 

For concision, we present only the verification of full-/ COBBLES, except for con- 
servation properties, for which it is revealing to compare full-/ and 5f approaches, 
and for verification of drag and diffusion collision operator, since reference material is 
based on hypothesis of 5f model. As a preliminary test, we compared linear growth 
rate and real frequency measured in COBBLES simulations with those obtained with 
Davies solver in the Krook case, or by solving the eigenvalue problem in the Fokker- 
Planck case, and found good quantitative agreement, which is illustrated in Table 
[331 



3.4.1 Collisionless closed system (7^ = t'a./.d = 0) 

Our purpose is to test nonlinear capabilities of COBBLES. Let us consider the simpler 
collisionless Vlasov-Poisson model without external damping, corresponding to the 
BB model without any collision nor extrinsic dissipation. In the unstable case, linear 
growth goes on until a significant number of resonant particle trajectories are modified 
by electrostatic trapping. In the nonlinear phase, the distribution develops an island 
structure in phase-space, and becomes flat on average in the resonant velocity region. 
The instability saturates and linear theory breaks down. As a measure of the electric 
field amplitude Eq, we use uj^, the bounce frequency of particles that are deeply 
trapped in the electrostatic potential, which is defined by uj^ = kEo. O'Neil extended 
the theory of collisionless wave-particle interaction in the nonlinear phase iO'N65j . 
within the assumptions "fL/'-^b ^ 1 and uj/ujh 3> 1. He obtained an analytic estimation 
of the evolution of wave amplitude. In the small-time limit, uji,t <^ 1, the electric field 
amplitude is estimated as 



uJb{t) 7l f 2ujbt\ 
— -— = exp / dK [ 1 — cos . 



(3.64) 



Fig. 3.7 shows the evolution of normalized bounce frequency uJb/^L, along with 
snapshots of the distribution function, for initial distribution B. We recover the linear 
growth rate obtained from Davies solver within 1% error. The nonlinear evolution of 



the wave is in qualitative agreement with the analytic estimation (3.64) in its validity 
limit (for the first few amplitude oscillations). Although it is impossible to quanti- 
tatively compare all the features of this analytic solution because of an ambiguity in 
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Fig. 3.7: Full-/ COBBLES simulation with — Va.f,d — 0, for initial distribution B, 
with X Ny = 256 x 2048 grid points, (a) Nonlinear evolution of normalized bounce 
frequency, (b) Snapshots of distribution function. 



the initial time in Eq.(3.64), we observe a good agreement for the amplitude oscilla- 
tions frequency, and for the relative amplitude of these oscillations compared to the 
saturation level. Furthermore, the saturation level is close to the value uJb/^L ^ 3.2, 
which was numerically obtained in Refs. |CD93| lOWMTlj with the 5f BB model. 

In the time- asymptotic limit, assuming some infinitesimal amount of collision, 
the steady-state solution of the Vlasov-Maxwell system is a distribution given as a 
function of the energy only. This BGK solution is consistent with a non-zero electric 
field. Fig. |3.8| is a contour plot of the distribution function in the time-asymptotic limit 
of numerical simulation, on which several constant energy curves are superposed. We 
clearly observe an island structure, which agrees with the BGK solution. This island 
is topologically different from the initial condition, thus some collisions are needed to 
violate Liouville's theorem and obtain the BGK solution. In numerical simulations, 
finite numerical dissipation, which is due to interpolation on a discretized grid, smears 
out fine-scale structures near the separatrix, allowing a reconnection of contour lines 
of/. 



3.4.2 Conserved quantities 

Total particle number in simulations is calculated by taking the sum over the compu- 
tation domain of the integrated value of the electronic distribution, N{tn) = j P?j- 
When Va is a constant, the relative error in particle conservation is, as expected from 
a locally conservative scheme, of the order of machine precision (We are working 
with 64 bits double-precision variables, which use 8 bits for the exponent and 56 bits 
for the precision, so that 2^^^ ^ 10^^ is the minimum numerical error). Even when 



Vaiv) has the velocity dependence of the equation (3.6), the relative error is negligible 
(< 10-^ %), as shown in Fig. Is^gta). 



In both cases, numerical simulations show good fidelity to the power balance, 
even for a relatively small number of grid points. The relative error in power balance, 



\Pe + Pd + Ph \ /{\Pe\ + \Pd\ + \Ph\), is included in Fig. [3^ A direct comparison 
between 5f and full-/ is not really meaningful, since simulation parameters, and 
definitions of Pe and Pd are different. Fig. |3.10| illustrates how the different power 
transfers (normalized to Pq = irvu^j^/k'^) compensate with each others. 

Entropy can be seen as a measure of numerical dissipation, since it grows as small 
structures dissipate. In full-/ simulations we define entropy as a sum over grid points 
of /log/, and we check that / is always strictly positive. In df simulations, we 
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Fig. 3.8: Contour plot of the time-asymptotic distribution function (solid curves) and 
constant energy curves (dashed curves) for the simulation of Fig. 3.7 



10'' 
10" 



5 10-° 

)_ 

05 10-8 
05 

1 10-^° 

cu 

°= 10-^2 



10 

10" 



■14 



Full-f'COB^LES ■ 


1 


'S- \ Error in — - 




^.1 total density 




\ Error in 




power balance 




Error in 




total entropy 


(a) 



20 



40 



5f COBBLES 



Error in 

total density 

Error in 

power balance 

Error in 

total entropy 

u hii jliL.j.i ■ iL , nL Ihf JId L ilk Jiiii I |»i 1 1 alii ■ I 



(b) 



60 



20 



40 



60 



Fig. 3.9: Conservation properties of COBBLES. Time-evolution of the relative error 
in total particle number calculated with from /, in total particle number calculated 
from p ("total density"), in power balance, and in total entropy, for steady-state 
simulations, (a) Full-/ simulation with initial distribution B, Va{v > u^) = 0.002, 
7d = 0.002, iVa; X 7V„ = 64 X 512 grids, and CFL = 0.9. (b) 5/ simulation with 
= 0.1, 7d = 0.05, Va = 0.05, iV^, x iV^ = 64 x 512 grids, and CFL = 1.1. 
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Fig. 3.10: Time-evolution of the power balance. The power transfer is normalized to 
Pg = 7rwfl7|,/fc^. Pe-, Pd and P^ are the electric field, dissipative, and particle power 
transfers, respectively. The solid curve is the sum of these three power transfers, (a) 



Power balance in a full-/ simulation with same parameters as Fig. 3.9 a), (b) Power 
balance in a 5/ simulation with same parameters as Fig. |3.9|[b). 



arbitrarily assume a bulk distribution as = 21/^^;^^! in the resonant region, where 
is the minimum of /(f . The factor 2 is to ensure that / = -I- does not 
take negative values because of perturbations in near its minimum. The relative 
error in the total entropy is included in Fig. |3.9[ 



3.4.3 Benchmark 

We consider five kinds of behavior for the time-evolution of the instability in the 
Krook case, and produce the behavior bifurcation diagram in the (7^, Vg) space. 
These behaviors are illustrated in Chap. |4] (Fig. 4.1). The category is obtained by 



an analysis of the electric field energy density £(t) and of the spectrogram of electric 
field. A numerical solution is defined as 

1. Damped: if the asymptotic-time limit of £{€) is zero; 

2. Steady-state: if the asymptotic-time limit of £(i) is finite; 

3. Periodic: if for large enough t there is a period r for which E(t -f r) — > f (t); 

4. Chirping: if there is a spectral component whose frequency significantly shifts 
in time. 

5. Chaotic: if E(fy is bounded, but does not satisfy one of the previous conditions. 

The categories 1., 2., 3. and 5. are defined in the same way as Vann [ VDR+03j . and we 
added a new diagnosis for the characterization of chirping solutions. Each numerical 
solution is systematically categorized by an algorithm based on a decision tree which 
is based on the one developed by Vann. We describe this algorithm in Appendix [D| 
As a benchmark of both COBBLES code and our categorization algorithm, we re- 
produce results presented in Fig. 3. of Ref. jVDR+03] (Note that our definition of 7^ 
is consistent with Berk and Breizman's literature, and differs from Vann's article by a 
factor 2). The initial distribution is A, for which 7^ = 0.1981. The field energy of the 
initial perturbation is 2 x 10^* of the total energy, which corresponds to uJb/lL — 0.05 
at i = 0. We perform a series of simulations in the parameter space (7^, Va) , where 
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Fig. 3.11: Behavior bifurcation diagram. The classification of each solution is plotted 
in the (7^, Va) parameter space. The solid curve is the linear stability threshold 
obtained by solving the linear dispersion relation numerically. The parameters of 
these simulations are TV^; x iV^ = 64 x 512 and CFL — 1.2. 



Va is chosen as velocity-independent. We set the time-duration of each simulation to 
^max — 3000. In the categorization algorithm, we choose imin — 1000, ei = 10^^^, 
€2 = 0.05, £3 = 0.01, 64 = 10^^, and £5 — 0.25. The resulting behavior bifurcation 
diagram is shown in Fig. |3.11| The 1416 simulations used for this plot took approxi- 
mately 115 CPU hours on an Altix3700Bx2 array of Intel ltanium2 processors. The 
categorization of 92 % of these time-series is in agreement with the reference, most 
of the difference coming from a different way of sorting out chaotic from periodic 
solutions. This result is a further indication of the validity of both COBBLES and 
categorization algorithms. 



3.4.4 Steady-state with drag and diffusion 

To verify our implementation of collision operator with drag and diffusion, we con- 
firmed that a Gaussian perturbation in the velocity distribution follows the analytic 
solution of the diffusion equation in the absence of electric field and drag, and is sim- 
ply adverted at a rate v — v'j/k in the absence of electric field and diffusion. 

As an additional test, we compare nonlinear steady-state solutions between Sf- 



COBBLES and analytic predictions derived in Ref. [LilOO . Fig. 3.12 shows steady- 
states in (5/-C0BBLES simulations with different collision frequencies. Fig. 3.12[ b) 
in this manuscript, and Fig. 6.13 in Ref. |Lil09j . which share the same normalization, 
can be directly compared. We confirm quantitative agreement with theory. 

3.4.5 Multiple-modes interaction 

Though this feature is not used in this thesis, we also test multiple-waves capabilities 
of COBBLES. When many electrostatic waves are excited, the amplitude of each 
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Fig. 3.12: Steady-state saturation with drag and diffusion, df simulation with 7/^0 = 
0.02, 7d = 0.01, Vd = 0.1, and given in the legend, (a) Saturation level, (b) 
Saturated distribution function. 



wave grows exponentially until nonlinear saturation occurs, and each wave develops 
an island structure in phase-space. If the width of each island is much smaller than 
the distance between the phase velocities of two neighbouring waves, we can treat 
the problem as a superposition of the former single wave-particle problem. However, 
if island structures overlap each others, particle trajectories are not integrable. We 
consider a situation where there exists a velocity interval within which the phase 
velocities of many waves are close enough and their islands overlap. We perform a 
full-/ COBBLES simulation, without collisions nor external damping, with ub — 0.05, 

VTM = VTB = 4.0, VB = 16.0, t^max = -"min 30, L = 512, X = 512 X 64, 

initializing 10 waves with wave numbers km = 21:771/ L and a random phase for each 



wave m. Fig. 3.13 shows the position and width of each island, and trajectories 
in the velocity direction of three test particles evolving within the resonant region. 
We observe overlapping of islands, and resonant particles seem to undergo Brownian 
motion in the velocity direction. 

When particle diffusion time is much longer than correlation time, quasi-linear 
theory |SG69| predicts velocity diffusion of particles in the resonant region, leading 
to a flattening of the distribution, as we observe in numerical simulations. In the 
resonant region, the quasi-linear diffusion coefficient -Dql can be estimated as 

i?QL - , (3.65) 

where E}^^ is the Fourier component for the wave number km of the electric field, and 
Au/j is the width of the whole resonant region. 

Another way of estimating the diffusion coefficient involves the variance of the 
displacement in velocity of a large number of test particles. For any time interval At 
larger than the correlation time, but smaller than the distribution relaxation time, 
this estimated coefficient Z3p is given by 

{[v{ta + M) - vito)]") 
2At 

where angle brackets represent an ensemble average. 

In our simulation, we estimate Dp by following the trajectories of 3 x 10^ test 
particles, which are initialized with an uniform distribution over the resonant region. 



and a random position. Fig. 3.14 shows that -Dql matches Dp in simulation, even as 
we double At. 
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Fig. 3.13: Evolution and overlapping of multiple islands. Dashed and dotted lines are 
trajectories of three test-particles. 




40 



Tab. 3.4: Analogies between the BB model and the reduced model for the TAE. 



Harmonic oscillator 



BB model 



ID model for TAEs 



Hamiltonian 
^DP + V cosV' 

Action 
/ 

Angle 

Effective mass 
D 

Oscillations amplitude 
V 



Hamiltonian 

1{v-vk)^ + ^(x,t) 

Velocity in the wave frame 
I{v) = (w - V'R)/k 

Position in the wave frame 
il>{x, t) ~ kx — Lot 

D = P 

Electric field amplitude 
V = u:l/k^ 



Hamiltonian 

\DP + V cosi/' 

Deviation from resonant surface 



+ nQ — ujt 



D : 



'Ho 



Magnetic field amplitude 
V = ~ieC{t)VpiJgR, Jcr) - SA_ 



3.5 The BB model as a paradigm for the TAE 

The BB problem can be put in Hamiltonian form with the Hamiltonian given in 
Eq. ( 3.13| ). Let us make a canonical transformation with the generating function 



S = VB.ix — v-Rt/2) to a moving-frame coordinate set, ('0, I), where ijj = kx — Lot and 
I = [v — VYi)/k. The new Hamiltonian, 

^2 ^2 ^2 

= h - Ilo - ^ ^ —P + ^ cosV, (3.67) 
takes a standard form, which is shared with the effective Hamiltonian of the TAE, 



Eq. (2.69). Therefore, the BB problem is a simple ID model that is homothetic to a 
whole class of instabilities, including EP-driven TAEs. 

Fig. |3.15| is a schematic representation of wave-particle interactions relevant to 
TAEs. A dotted-dashed rectangle represents the physics encompassed by the BB rep- 
resentation of TAEs. All physics outside this rectangle are treated as input parameters 
in the BB model. 

Table [3^ summarizes the parallels between the BB model and the reduced model 



for the TAE. 



41 



3 




m 



—\ 

3 
o 

Q. 
fD 

in 



^Kinetic ^ 



effects 




o 



Particule 




source 
Particule 



Loss 



sink 
Collisions 



Q. 

3 

■D 

in 



Q. 
C 





H 









fD 



















3 



fD 



fD 


H 


fD 




n 


fD 
—\ 










3 






C/1 





orbit 



(/I 3! 

rT I 
fD Q) 



Slowing-down 



n2 

- o 
00 n 

fD 



Fig. 3.15: Schematic representation of wave-particle interactions relevant to the TAE. 
The dotted-dashed rectangle represents the limit of the BB representation of TAEs. 
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Chapter 4 

Kinetic nonlinearities 



To different nonlinear behaviors of the BB model correspond different regimes of res- 
onant energy transfer. In terms of alpha-particle issues, there is an interplay between 
two effects. On the one hand, energy transfer from particles to the wave, which is sup- 
ported by the thermal plasma, is favorable since the energy from fusion reactions must 
be channeled to heat the bulk plasma. On the other hand, a large wave amplitude 
is unfavorable, since it is associated with transport and energetic-particle ejection. A 
survey of kinetic nonlinearities provides important insight into the optimum regime. 



In Sec. 4.1 we perform a systematic parameter scan in (7^, Va) in the full-/, Krook 
case. We confirm the validity of available theory, and show limits of Sf approach. In 
Sec. |4.2| we investigate nonlinear chirping features, since they can provide precious 
information about the state of the plasma. Existing quantitative predictions of these 
features are verified for both collision operators, and we extend theory by including 
the effects of beam distribution shape, finite collision frequency, and drag. In Sec. |4.3[ 
we investigate instabilities that arise in a regime where linear theory predicts wave 
damping, provided that initial perturbation is large enough. We propose a mechanism 
to explain the apparent contradiction between linear theory and the behavior of sub- 
critical instabilities observed in simulations, and perform a numerical investigation of 
an initial amplitude threshold. 



4.1 Nonlinear regimes 

Theories [BB90( IBBP96t IBBP97b| IBBP98| lBBC+99| have been developed by Berk, 
Breizman, and coworkers, to quantitavely predict nonlinear behaviors in various pa- 
rameter regimes and to explain underlying mechanisms. On the one hand, some of 
these theories have been validated by numerical simulations based on the 6f model. 
A concern with this perturbative approach is that, as instability grows, the resonant 
region may expand and ultimately include a significant portion of bulk particles. An 
other concern is that, when chirping occurs, corresponding resonant velocity may 
propagate into the bulk. In such situations, kinetic effects of bulk plasma should also 
be taken into account. On the other hand, full-/ simulations have been performed by 
Vann and coworkers [VDR^OS] . However, as we show in 3.3.2 this approach shows 



some difficulty in simulating situations considered in the aforementioned theories, 
which assume a plasma near marginal stability with a cold bulk and a weak beam. In 
fact, these theories have not been validated with this approach. Filling the gap be- 
tween these two fronts of the current state of research, with quantitative comparisons 
between available theory and full-/ model, is the aim of this section. We investigate 
the validity of analytic theories for the following nonlinear features, 

• saturation level in a parameter regime above marginal stability; 
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• saturation level and bifurcation criterion between steady-state and periodic so- 
lutions near marginal stability; 



time-evolution of a frequency shifting mode. 



We choose initial bunip-on-tail distribution B, which was introduced in 3.1.2 The 
parameters of distribution B were actually chosen so that we stay within the validity 
limit of these theories, with sufficiently cold bulk and sufficiently weak warm beam. 
We recall these parameters as ub = 0.1, vtp — 0.2, vtb = 3.0, vb = 5.0, which give, 
for k = 0.3, jL = 0.0324 and uj = 0.925. As mentioned before, i^a is a function of the 
velocity such that collisions affect only the beam particles. The field energy of initial 
perturbation is 2 x 10^^ of total energy, or i^b/lL = 0.3. 



4.1.1 Nonlinear saturation 

Fig. |4.1| shows four examples of nonlinear saturation in the unstable case, correspond- 
ing to steady-state, periodic, chaotic, and chirping behaviors, obtained by varying 
Va at fixed 7^. Spectrums are obtained by applying Fast Fourier Transform to time- 
series, which are filtered by a Hann window [PT VF92. . Included are both spectrum of 
electric field amplitude, where we consider only times after nonlinear saturation, and 
spectrogram of the electric field measured at some arbitrary point in configuration 
space. 



4.1.2 Scan in the (7^, Ua) space 

For the benchmark in |3.4.3[ we set a same value for the maximum time of every 
numerical simulations. However, we must now take into account computational cost, 
which is much larger because of a decrease of by one order of magnitude. As 
we approach marginal stability, the time window must be increasingly large to suc- 
cessfully capture the nonlinear behavior. To reduce computational cost, we choose a 
time-window size as a function of 7, as 

27r 

tmax = 20^. (4.1) 

l7l 

The frequency of amplitude oscillations is of the order of ujb, which is empirically 
of the order of 7 after the transient phase, so that such time windows contain at 
least a few amplitude oscillations, enough to sort steady-state, periodic and chaotic 
responses. In the categorization algorithm described in App. [D] we choose i„ii„ = 
^max/2, and each time series is sampled every Atg — 20. cq = 10^^** is chosen as a 
free-streaming criterion, eg = 0.05, ey = 0.05 are used to sort out chirping solutions, 
and the other e-thresholds are the same as above (ei = 10"-^^, £2 = 0.05, £3 = 0.01, 
£4 = 10~^, and £5 = 0.25). The behavior of wave amplitude time-series obtained by 
full-/ COBBLES is characterized in Fig. |4.2[ The 391 simulations used for this plot 
required approximately 15000 CPU hours on an Altix3700Bx2 array of Intel Itanium2 
processors. 

Agreement between linear stability threshold and the boundary between linearly 
stable and unstable simulations confirms that the problem of recurrence is taken care 
of by the free streaming test in our categorization algorithm. When <^ 7^ and 
Va ^ 7cij a bursty behavior, characterized by a succession of bursts with characteristic 
growth and decay rates of and 7^, respectively, and with a quiescent phase in 
between that lasts a time l/va, as described in |BBY92j . is expected. A few solutions 
in the chaotic region appear to follow this picture. However, most of the chaotic 
solutions do not feature a significantly quiescent phase. Consequently, an attempt 
at sorting out a pulsating regime from the chaotic region seems vain. For small 
collision rates, we observe instabilities in the linearly stable region, which suggests 
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50 100 
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Fig. 4.1: Typical nonlinear behavior in steady-state (a-c), periodic (d-f), chaotic (g-i) 
and chirping (j-1) cases. Full-/ simulations with initial distribution B, ^4 = 0.03, and 
values for Va are 0.02, 0.008, 0.005, and 0.00002, respectively. Left: Time-evolution 
of electric field amplitude. Center: Fourier spectrum of electric field amplitude after 
t = 3000, with a time- window At = 10^. Right: Spectrogram of electric field at a; = 0, 
with a moving Fourier-window of width At = 400. 
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Fig. 4.2: Behavior bifurcation diagram for a cold bulk, weak warm beam distribution. 
The classification of each solution is plotted in the (7^;, i>a) parameter space. The 
solid curve is the linear stability threshold obtained by Davies solver. The parameters 
of these full-/ simulations are CFL = 3.0, Umin = —10, Wmax — 18, Ny = 2048 and 
Nx ranges from 128 to 256. Smaller diamonds and triangles on the right of the linear 
stability threshold represent subcritical instabilities. 
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Fig. 4.3: Saturation level at a given collision frequency VallL = 0.1 for distribution B 
(L0W-7L case), and for higher beam density and lower beam temperature (Higher-7i 
case). The parameters of full-/ simulations are those of Fig. 4.2 
simulations are such that 7^ take the same two values. 



Parameters of (5/ 



the possibility of subcritical instabilities. This effect is discussed in Sec. |4.3| The 



chirping regime is discussed in details in Sec. 4.2 The physics of several other regions 



of this diagram is discussed in the remaining of this section. 

4.1.3 Steady-state above marginal stability (7,^ ~ ^ 7l) 

When external damping and distribution relaxation are of the same order and both 
are small compared to the linear drive, we expect and observe the saturation of wave 
amplitude to a steady-state in the time-asymptotic limit. To estimate a saturation 
level, we assume a rate of annihilation of beam particles much smaller than the sat- 
urated bounce frequency, i^a ^ Wfc at i — > cxo. We also assume that the resonant 
region is narrow compared to the resonant velocity, 4wft/fc <C w/Zc, so that we can 
assume that the contribution to resonant power transfer comes from a narrow region 
around vr- Berk and Breizman derived a relation yielding the saturation level in this 
situation |BB90| . 

uji, = 1.96 — 7i. (4.2) 

Id 

Thus, if we re-normalize all quantities to the linear growth rate, then within the 
aforementioned assumptions, the saturation level depends only on the ratio of Va to 
Id- 

We investigate the validity of this theory by numerically computing the scaling 
law for the saturation level at a given normalized relaxation rate Va/lL — 0.1. In a 
previous work [BBPQSaj . such a scan has been done using a 5f particle code, and the 
results showed good agreement with analytic prediction in a region where 7^ ^ Va- 



Fig. 4.3 shows the saturation level obtained from theory, Eq.(4.2), and from both 5f 
and full-/ COBBLES simulations. When the initial distribution is distribution B, we 
observe quantitative agreement between theory and both 5f and full-/ simulations in 
the parameter region jd ^ i^a- 

To reveal some limitations of Sf model, the same computation is done for a dis- 
tribution with a slightly higher beam density, ns = 0.15, and a slightly lower beam 
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temperature, vb = 2.5, giving 7^ = 0.067 instead of 0.032. When we use Sf model, 
the scaling law is roughly independent of and is in agreement with theory in a 
parameter region jd ~ i^a, in agreement with aforementioned work. On the other 
hand, when we take into account the evolution of bulk plasma, we observe a signif- 
icant dependency of the saturation amplitude on the linear growth rate. For larger 
7i, we find some discrepancy with theory in the low region, because the island 
width Aw becomes of the order of the resonant velocity {Av/vr = 0.14a;b/7i, in the 
I0W-7L case, and Av/vr = 0.29ujb/jL in the higher-7L case). This result shows that 
it is necessary to take into account the effect of bulk particles to accurately discuss 
the validity limit of this theory. 



4.1.4 Near-marginal steady-state and periodicity (7 ~ 7l — 

Id < 7l) 

When 7 ^ 7l, a reduced integral equation for the time evolution of electric field 
amplitude has been developed using an extension based on the closeness to marginal 
stability [BBP96j. Within the assumption Wfc/7 <C 1, 

"I ^ Jt/2 Jt-ti 

For a cold bulk, warm beam distribution, in the collisionless limit, as we approach 
marginal stability, Eq. |3.44| reduces to 

7 ~ 7lo - Id, (4.4) 



which agrees with the linear part of the latter integral equation (4.3 1. In Ref. |BBP96) 



the analytic treatment is carried on by normalizing time by jlo ~ Id- 



We observe that the relation (4.4) is a good approximation in most of the param- 
eter space. However, as we get closer to the linear stability threshold, the relative 
error |7Lo~7d~7|/(|7Lo~7d| + |7|) approaches unity for finite collisions. In addition, 
for our choice of distribution, there is a 14% discrepancy of 7^0 = 0.0368 compared 
to 7l = 0.0324. We infer that we can replace 71,0 ~ 7d by 7 in the integral equation 



(4.3) and use 7 itself as the relevant choice of normalization parameter. 
This procedure yields a steady solution. 



L^i = 2^/2^.2 /JL. (4.5) 

V 7-LO 

A series of simulations near marginal stability (0.005 < 7/7L < 0.02), for fa spanning 
2 orders of magnitude, confirms the validity of the latter expression. Fig. |4.4| shows 
quantitative agreement with the saturation level of numerical solutions. 



Nonlinear stability analysis reveals that the steady solution (4.5 ) is unstable when 



>^a < i^cr, with j/cr = 4.47. To assess this criterion for the bifurcation from steady- 



state to periodic solutions, a zoom in the behavior bifurcation diagram (Fig. 4.2) in 
a region near marginal stability where this bifurcation occurs is presented in Fig. 4.5 
We observe a qualitative agreement between the steady-periodic boundary and I'cr- 
However, when 7/7^ < 0.01, chaotic solutions appear for ^ Vcr- This discrepancy 
is explained by the existence of nonlinear excitations. As we approach marginal sta- 
bility, the nonlinear behavior becomes sensible to the initial perturbation level. To 
prove this point, we perform a series of simulations in the vicinity of the bifurcation 
with an initial amplitude reduced from ujb/lL = 0.3 to ijJb/^L = 3 x lO"''. Fig. 4.6 



shows the values of v/^ for the bifurcation between steady-state and periodic solu- 
tions. The bifurcation occurs somewhere in between. We confirm that fcr/7 stays 
close to the predicted value of 4.4 for smaller values of 7. 
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Fig. 4.5: Zoom in the behavior bifurcation diagram (Fig. 4.2 ) in a region near marginal 
stabihty where steady-periodic bifurcation occurs. Dashed Hue is the critical distri- 
bution relaxation. Smaller diamonds on the right of the linear stability threshold 
correspond to subcritical instabilities. 
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Fig. 4.6: Critical distribution relaxation for steady - periodic bifurcation near marginal 



stability obtained from full-/ simulations with numerical parameters of Fig. 4.2 



4.2 Nonlinear features of chirping 



In the collisionless limit, when J^a < 7 ^ 7l, the integral equation (4.3) is consis- 
tent with explosive solutions that diverge in a finite time, which suggests that the 
mode energy is partitioned into several spectral components. The resulting sideband 
frequencies have been observed to shift both upwardly and downwardly [BBP+QTc] . 
the frequency shift (5w(t) increasing in time. 

These chirping solutions arise when hole and clump structures |BBP97a] are 
formed in phase-space. They belong to a chaotic regime, and each chirping event 
is slightly different. In this section, we are interested in the nonlinear chirping char- 
acteristics, averaged over a significant number of chirping events. In particular, in 
our simulations, the first chirping event is observed to stand out from the statistics, 
with a larger extent of chirping - up to twice as much as any other one of the follow- 
ing series of repetitive chirping. This may be due to the fact that the first chirping 
benefits from a perfectly constant velocity-slope, while following events suffer from 
the interference of phase-space structures that remain from previous chirping events. 
Since the latter condition seems more experimentally-relevant, the first chirping is 
ignored in the present analysis, unless stated otherwise. 

Since we want to use chirping features as experimental diagnostics, it is necessary 
to validate and develop corresponding theory. These features are quantified from raw 
simulation results in this section, and from experimental data in Chap. [Sj using an 
algorithm described in Appendix [E] 

4.2.1 Holes and clumps 

Holes and clumps are nonlinear coherent structures with time-dependant velocities. 
In general, several holes and several clumps, with different amplitude, coexist, as 
shown in Fig. |4.7| |a), which is a snapshot of the velocity distribution for a Krook 5f 
simulation. In Fig. |4.7[ b) , we plot D as a function of the real frequency w and the 
growth rate 7, where 

Z?(u;,7) ^ 7 + 7d - - ^ / ^ rdt', (4.6) 



2fc Jr (7 + i^a) + « (fc w - u) 



so that _D = gives the linear dispersion relation Eq. (3.58). This plot suggests that 



to each hole/clump corresponds an eigenmode, with frequency and growth rate both 
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Fig. 4.7: Hole and clumps, (a) Snapshot of velocity distribution for a Sf simulation 
with = 0.1, 7<i = 0.05, Va = 0.001, iV^ x iV^, = 64 x 2048, and CFL = 1.1. 
(b) Determinant of the linear dispersion relation, where the distribution function is 
taken from the simulation a,t t — 600. Logarithmic color scale spanning 6 orders of 
magnitude. The center of each dark area corresponds to a linear mode. 
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functions of time. 



4.2.2 Chirping velocity (or sweeping rate) 
Available theory 

Ref. |BBP97bj shows how one can isolate one spectral component and model it by a 
BGK wave to obtain the time-evolution of one chirping event. This theory is based 
on the following assumptions: 

• The resonant velocity of a hole/clump evolves slowly enough for trapped particle 
orbits to keep their coherency, Slo/ujI^ 5uj/ujI ^ 1; 

• The width of a hole/clump evolves slowly enough for trapped particle orbits to 
keep their coherency, Wft/o;^ ^ 1; 

• Holes and clumps are narrow enough that they don't overlay each others, uJi,/Sio ^ 
1. 

Within the above assumptions, the perturbation of passing particle distribution is 
negligible, and a bounce- average treatment of trapped particle distribution yields the 
frequency shift, in the collisionless limit, as 

Soj{t) = a-fLoVjdt, (4.7) 

with a w 0.44 ; and a saturation level as 

ujb ~ 0.54 7io. (4.8) 

These analytic expressions have been found to agree with ID Sf particle simula- 
tions, pBP97bj . with both Krook and diffusion-only collision operators, and with 3D 
HAGIS simulations |PBG+04| . 

Numerical validation 

When is finite and Va is small enough, we observe such chirping solutions in full-/ 
COBBLES simulations. Fig. |4.8| shows the time evolution of field amplitude, with 
initial distribution B, when 7^ = 0.035, and Va — W^*, so that 7 — 0.05 7lo- The 



simulation result agrees with a numerical solution of reduced equation (4.3), until the 
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field amplitude approaches the applicability limit, 
close to analytic prediction Eq. (4.8). 
Fig. 



After saturation, the solution is 



4.9 



a) shows the spectrogram of electric field. Frequency sweeping events 
occur repetitively as the slope of the distribution is successively recovered and flat- 
tened. The time evolution of each chirping event is slightly variable. To quantify the 



agreement with theory, we extract the 12 largest upshifting branches. In Fig. |4.10 



we show these branches shifted to the same initial time, where each point is obtained 
by interpolating local maxima from the discrete frequency spectrum. From this plot, 
we conclude that, before border effects occur, chirping evolution follows a square-root 
law in time, as expected from theory. This suggests a possibility to recover the prod- 
uct 7lo \/ld from such power spectrum. By fitting a square-root law to the branches 
of Fig. 4.10 we obtain an average value of 7lo >/7d = 0.0057, with a standard de- 

0.0061. When a plasma is close 

' ~ 7d^^)j and this chirping 
diagnostics may be used as a rough estimation of the extrinsic damping rate of the 
bulk plasma 7^. The agreement of the expression Eq. (4.7) with chirping found in 
experimental devices [PBG"'"04] supports this claim. 

As we consider a single mode with a fixed wave number, sweeping frequencies 



) we obtain an average value of 7^0 \ 
viation of 14 %, when the input value is 7^0 \fld 
to marginal stability, we can assume 7^0 ~ 7d (7lo 



correspond to evolving structures in velocity distribution. In Fig. |4.9[ b), we observe 
formation and evolution of hole/clump pairs in phase-space, and they show clear 
correlation with peaks of the spectrogram shown in Fig. 4.9 |a). It should be noted 



that near marginal stability, 7^ > 0.47^ is given as a necessary condition for hole- 
clump pair creation in Ref. BBP+97C] . Although there is yet no theory in the opposite 
limit 7 7/,, in our simulations we observe frequency sweeping even when 
A spectrogram is shown in Fig. 



4.11 



for 7d/7L = 0.2, Va/lL = 0.02. Althought 
the chirping is not as pronounced as in Fig. |4.9[ a), we observe that the dominant 
frequency sweeps 5% of its initial value. 



Non-adiabatic chirping 

In Chap, [sj we consider a regime with relatively fast sweeping, 5uj/uj^ ss 0.5, which 
approaches the limit of validity of the above theory. 6lo juy^ can be seen as a measure 
of hole/clump adiabaticity, and is roughly proportional to (^^dMa)^^'^ llm in the Krook 
case. When Slj/ljI sa 0.5, 4wfc/(5a; « 27r/ajb, in other words a hole or a clump is shifted 
by its width in a bounce time of deeply trapped particles. In this regime, the previous 
analytic treatment is not relevant. 

Clarifying this point is easier if we avoid effects of bulk evolution and of the shape 
of beam distribution. Thus we switch to the 5f model. Fig. |4. 12] shows spectrograms 
of chirping 5f COBBLES simulations, first in a regime which satisfies the assumptions 
of the above theory, then in a regime out of its scope. We observe similar square-root 
dependency of the frequency shift in time. This suggests that we can introduce the 
effect of non-adiabaticity on chirping velocity as a correction parameter fi, defined as 



(4.9) 



which is a priori a function of all input parameters 
iwl^ = 0.1 in Fig 



4.13 



P is obtained numerically for 
Results show that chirping velocity slows down compared 
to theory as we leave the adiabatic limit. We confirm that, inside the validity limit 
of the above theory, /3 approaches unity. Even for relatively large values of Suj/uJ^, 
chirping velocity has a smooth dependency on the kinetic parameters. The latter 
point is crucial with regard to the validity of the procedure described in Chap. [5j 
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500 1000 1500 

Fig. 4.9: (a) Logarithmic scale plot of the time evolution of the frequency power 
spectrum. At each time, the spectrum is normalized to its maximum value, (b) 
Logarithmic scale plot of the time evolution of 5/ = / — /q. 
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Numerical result 

Analytic theory 
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Fig. 4.10: Time-evolution of the 12 largest upshifting branches in logarithmic scale. 



Dotted line is the first chirping event. Dashed line is the analytic prediction Eq. (4.7 1 
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Fig. 4.11: Logarithmic scale plot of the time evolution of the frequency power spec- 
trum for 7d = 0.2 7l, Va = 0.02 7^. 
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Fig. 4.12: Spectrogram of the electric field in Sf simulations, obtained with a moving 
Fourier window of size A = 510, for jlo = 0.1 and (a) jd/lLo — 0.4, ^d^a/ltf) — 0.008, 
/3 = 1.0; (b) 7d/7L0 = 0.9, ^dVa/llo = 0.043, P = 0.57. Solid lines show the chirping 
velocity predicted by Eq. (4.7), with correction coefficient /3, and chirping lifetime 



predicted by Eq. (4.22 1. In the remaining of this thesis, the logarithmic color scale 



for each spectrogram spans 3 orders of magnitude. 
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Fig. 4.13: Correction to Eq. ( |4.7[ ) when the time-scale of frequency shift is relatively 
short compared to the bounce period, measured in 5f simulations with 7^0 = 0.1. 
The ratio ^d/^cs is such that 5uj{l/vcii) = 0.2 in Eq. (4.7 1, (a) In the Krook case, 
spectrograms for two extreme points are shown in Fig. 4.12 (b) With drag and 
diffusion, I'f/vd = 0.1. 
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Fig. 4.14: Effect of finite Krook collisions on chirping velocity. Spectrogram of Sf- 
COBBLES simulation with 7^0 = 0.05, 7<j = 0.045, = 4 x IQ-^^ x = 
128 X 4096, and CFL 1.3. Two straight lines correspond to Eq. ([ij]). Two bended 



lines correspond to Eq. (4.20). In both cases, the lines span a time interval 1.1/i^a 
[See Eq. (4.22)], and we have included a correction coefficient /3 — 1.23. 



Inset: 



on the beginning of the first chirping event. 



Effect of finite Krook collisions 



In Ref. |BBP97b] . Eq. ( |4.7[ ) is derived by changing variables to action-angle of the 
bounce-motion of particles trapped in a hole or a clump, and by noting that the 
unperturbed part /o(w) of / does not contribute in the resonant power transfer of 
the power balance, Eq. ( |3.30 1. Thus, the wave equation involves only the deviation 



from the unperturbed distribution at the center of the evolving hole/clump, g 
/ — /o(w + 6uj). Then, g is expanded in powers of e, 

g ^ go + egi + ... (4.10) 

where e = max((5aj/w^, Wb/w^, Wb/Jw) 1, and it is shown that to zeroth order in e, 
the real and imaginary parts of the wave equation are reduced to 

27L0 r^'— sj^ 

T^^ldufoii^o) Jo 



7. - -rrr^TT^ :Ti^odJ, (4.ii) 



Slu = 2?§V^ / (cosV)goclJ, (4.12) 



7ra;^9c^/o('^o) Jo 

where angle brackets indicates a bounce-average, ip is the spatial coordinate in a frame 
moving with the hole/clump, J is the bounce- motion action, Jmax = 8w(,/7r, and go 
is obtained by bounce-averaging the kinetic equation. Then, in the reference, the 
collisionless limit is consider, in which case go is simply 

50 (t) = foH - /o(c^ + M- (4-13) 



Eq. (4.7) follows by assuming a constant gradient for /q. 
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The above theory in the adiabatic hmit, and the above numerical investigation of 
(3 are vahd on a timescale smaher than a colhsion time. In the following, we include 
the effect of finite collisions, to explain a deviation from square-root law in longer 
timescales. In the Krook case, the bounce-averaged kinetic equation for g at lowest 
order is 

dgn , dfo 



dt 



dt 



whose solution is 



(4.14) 



go{t) = e-^"*/o(^) - foi^ + Sio) + Va / ^ 8i^{t')Y-^'' dt' . (4.15) 



As in the reference, we assume a constant gradient for /q, 

/o[c^ + <5c^(t)] - /o(c^) + Uoj)5uj[t) 

which yields 

50 (i) = ./o(^) 



(4.16) 



(4.17) 



''»(*'-*)(5a;(i')dt' - 5uj{t) 
We substitute the latter solution into Eqs. ( 4.1l|4.12 1 to find, in the limit dui/uil <^ 1, 

(4.18) 
(4.19) 



7d 
(5a; 



= 3 



,2 ' 



I67 



LO 



37r2 



1 - 



6uj{t') dt' 



We solve the latter equation system by expanding both ujh and Slu in powers of i^at- 
This lengthy but straightforward procedure yields 



6u!(t) = ±af3-fLQVldt 

2 



1 - -(Vat) + —(Vatf ^—(Vatf 

' ' 1890^ ' 



UJb{t) 



I67LO 

37r2 



1 



315 



(4.20) 
(4.21) 



where we have included the effect of non-adiabaticity. Note that Eqs. (4.7)-(4.^ 



are 



recovered in the coUisionless limit. The effect of finite collision is to reduce chirping 
extent by bending shifting branches. This effect is not negligible since 5oj is reduced 
by 27% after a collision time, which is of the order of chirping lifetime as we will see 
in |4.2.3[ This is illustrated in Fig. |4.14[ which is the spectrogram for a 5f Krook 
simulation. The beginning of the first chirping event is fitted to the coUisionless 
analytic prediction, Eq. (4.7), to determine the correction parameter /3, then we show 
the time-evolution of Eq. (4.20) with the same correction parameter, which is in good 
agreement with observed bended chirping. 



4.2.3 Chirping lifetime 

The resonant velocity of a hole (a clump) does not increase (decrease) indefinitely. 
We define the lifetime r of a chirping event as the time in which the corresponding 
power in the spectrogram decays below a fraction of the maximum amplitude that 
is reached during this chirping event. The maximum lifetime Tmax is the maximum 
reached by r during a time-series, ignoring the first chirping event and any minor 
event. It is reasonable to assume that the island structure is dissipated by collisional 
processes, in which case the maximum chirping lifetime should be of the form 

'^max — ; (4.22) 
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Fig. 4.15: Maximum lifetime of a hole or clump, far from marginal stability, 
Id/lLVi — 0.5, and near marginal stability, ^d/lLo — 0.9. (a) With a Krook colli- 
sion operator. Crosses correspond to 5 J simulations with initial distribution shown 
in Fig. |3.2[ a), while triangles correspond to full-/ simulations, where the initial dis- 
tribution is a bump-on-tail with a Gaussian beam. In both cases, 7^0 — 0.05. A solid 
line corresponds to Eq. (4.22). (b) 5f model with drag/diffusion collision operator 
(for a linear distribution only), for two different values of linear drive. A dashed line 
corresponds to Eq. ( |4.23| with La 0.16, a solid hne to Eq. ( |4.24[ ) with id = 10. The 
drag is chosen so that it does not significantly alter chirping lifetime, Vf/i'd = 0.1. 
An absence of point means that we do not observe repetitive chirping before the end 
of simulation [t = 100000). 
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in the Krook case, and 



(4.23) 



in the case with drag and diffusion, when Vf <^ Vd, where La and are constant 



parameters. In Fig. 4.15 we plot the maximum lifetime measured in (5/-C0BBLES 
simulations, where the ratio ^d/lw is chosen as 0.5 and 0.9, i.e. far from and close to 
marginal stability, respectively. A quantitative agreement is found with Eq. (4.22 1, 
with La — 1.1, for spanning 2 orders of magnitude. With the diffusive collision 



operator, chirping lifetime agrees with Eq. (4.23) only for low collisionality. For high 



collisionality, diffusion affects the width of a hole or clump during the first phase 
of their evolution, namely drive by free-energy extraction, which in turn affects the 
decay by diffusion. Since chirping observed in experiments has a lifetime of the order 
of T '-^ 500, we adopt a semi-empirical law obtained by a linear fit. 



id 



0.5 



(4.24) 



with id = 10. No repetitive chirping is found near marginal stability for 0.05 < 7lo i£ 
0.1, though longer computations may reveal this possibility. 

In Chap. [5] Eqs. (4.22) and (4.24) are used as diagnostics for the effective col- 
lision frequency, thus it is important that these results are not too sensitive to the 
shape of fast-particle distribution, which is not measured accurately enough in ex- 
periments. To investigate this point, we repeat the previous analysis (in the Krook 
case), this time with an initial bump-on-tail distribution with a Gaussian beam (with 
full-/ COBBLES) instead of a constant gradient, or linear, beam (as is imposed in 



Sf COBBLES). Fig. 4.15 'a) shows that the agreement is kept, even if the shape of 



the distribution has a significant effect on the extent of chirping as can be seen for 



example in Fig 4.9 'a 



4.2.4 Chirping asymmetry 

It should be noted that, in Fig. |4.14| for example, we observe a slight asymmetry 
between up- and down- shifting branches, which is in disagreement with theory, since 
the Sf BB problem with Krook collisions is symmetric around the resonant velocity 
Vfi. We infer that the observed asymmetry is a numerical effect due to a difference 
in round-off errors between v < vr and v > vn domains, which is amplified in 
such chaotic regime. In a long-time averaged point-of-view, chirping observed in 6f- 
COBBLES simulations should be symmetric, though. 

In this section, we consider significant chirping asymmetry, rather than the above 
negligible numerical effect. 



Effects of drag 

The drag term in the collision operator has a counter-intuitive effect on chirping 
asymmetry. Since the effect of drag on any phase-space structure is to advect it from 
large to small velocities, one could imagine that the presence of drag would make 
down-chirping dominant. The opposite is observed in our simulations. In Fig. |4.16[ 
we show a symmetric simulation obtained with negligible drag, and two asymmetric 
solutions with sign ifica nt drag. For Fig. 4.16[ c), we changed the sign in front of 



{v'j /k)dyf^ in Eq. (3.7). An explanation for this paradox is proposed in Ref. [LBS 10) 



4.2.5 Chirping quasi-period 

In some parameter regimes, chirping arising from the BB model with Krook collisions 
is quasi-periodic, although the phase between two major chirping events is generally 
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Fig. 4.16: Effect of drag on chirping asymmetry. Spectrograms of Sf COBBLES 
simulations, with input parameters 71,0 = 0.099, 7^ = 0.028, Vd = 9.2 x 10^"^. (a) 
Natural drag, h'f = 5.4 x 10^"^. (b) Negligible drag, lyj — 8.0 x 10^^. (c) Inverted 
drag, Vf = 5.4 x 10^^. 




Fig. 4.17: (a) Time evolution of the bounce frequency. 6f simulations with 7^0 — 0.1, 
7rf = 0.05, and i/cS = 0.002. In the drag/diffusion case, jv^ = 0.3 (b) Distribution 
function, normalized to ^-iQl\f2T:. Solid line is the initial condition, dashed lines are 
for i/efft = 2. 
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Fig. 4.18: Periodic chirping regime. Spectrograms of 6f COBBLES simulations, 
with input parameters jlq = 0.099, 7^ = 0.028, and Vf/vj, — 0.2. (a) Less colhsions, 
Vf = 2.9x10-3, j/rf = 1.5x10-2. (b) More colhsions, Vf = 4.4x IQ-^, Vd = 2.2x10-^. 
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Fig. 4.19: Dragging regime. Spectrograms of Sf COBBLES simulations, with input 
parameters 7^0 = 0.099, 7^ = 0.028, and Vjjvd = 0.6. (a) Less collisions, Vj = 
9.7 X 10-3, = 1.7 X 10-2. (b) More collisions, Vj = 1.5 x lO^^, = 2.6 x 10"^. 



not quiet, but filled with minor chirping events. In a regime where Vf <C v^^ chirping 
arising from the BB model with drag/diffusion collisions is quasi-periodic too, but 



this time with clear quiescent phases in-between chirping events. In Fig 4.17 



which shows periodic decay and recovery of perturbation amplitude, corresponding 
to major chirping events, we observe qualitatively different behavior between the 
two collision models. The velocity distribution after nonlinear saturation shown in 
Fig 4.17| ^b) illustrates the fact that several holes and clumps with different amplitudes 
co-exist in the Krook case, while diffusion smooths out fine-scale structures, which 
explains isolated chirping events we observe in the drag-diffusion case (See for example 
Fig. |5.6[ c)). In both cases, no analytic theory has been developed that predicts the 
average time between two chirping events, Aichirp- However, conceptually, there exists 
some relation with a subset of the input parameters. This concept is useful to develop 
diagnostic for TAEs in Chap. [5] 

In the case with drag and diffusion, we found essentially two regimes, depending 
on the ratio Uf/vd, with a threshold r/^. When Vf /va < ffd, the chirping instability 
is in a periodic chirping regime as described above, with chirping events repeating at 
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regular intervals and a clear quiescent phase between two successive chirping events. 
This behavior is illustrated in Fig. 4.18 which shows two spectrograms with Vf /vd = 
0.2. In this regime, Aichirp decreases as drag is enhanced, as can be seen by comparing 
Fig. 4.18[ a) with Fig. 4.18 'b), since it accelerates the process of velocity-distribution 
recovery. However, Vf is not the only factor that determines Aichirp- If linear drive 
is decreased, more time is required to recover a gradient steep enough to overcome 
damping, hence Aicinrp increases. Diffusion tends to recover the initial distribution, 
but much less efficiently than friction. It also tends to counteract formation of holes 
and clumps, hence Aichirp slightly increases with increasing Vd- On the contrary, 
when Vf/vd > Tfd, dynamical friction significantly modifies the shape of a chirping 
branch, up to a situation where the direction of a chirping hole is reversed. This is 
illustrated in Fig. 4.19 which shows two spectrograms with Vj jvd — 0.6. We refer to 
this latter regime as dragging regime. Chirping such as in Fig. 4.19 ^a) has recently 
been observed and explained in Ref. 



[LBSI O . This is called as hooked chirping. 
In Fig. 4.19[ b), drag and frequency sweeping contradicting effects seem to balance, 
creating a hole with a stable frequency shift. Though we found a threshold ry^ ~ 1 
in our range of investigation, more systematic scan remains to be done. 

Note that in a case with diffusion only [va ~ I'f ^ 0), after a first chirping event, 
the perturbation is damped to zero, and the gradient of at resonant velocity is 
significantly reduced from initial condition. The original gradient does not recover in 
a reasonable time for typical experimental values of Vd- Drag is an essential ingredient 
of repetitive chirping, since, by adverting f — fo from positive to negative velocities, 
it replaces, at the resonant velocity, a plateau by the large gradient which was formed 
at slightly larger velocity compared to the plateau. 



4.3 Subcritical instabilities (7 < 0) 

Instabilities in a regime where the wave is linearly stable have been observed in Sf 
particle simulations [BBC+99) . We recover such instabilities with both df and full-/ 
COBBLES. Fig. 4.20| ^a) and (c) show the time-evolution of electric field amplitude 



obtained for distribution B, with jd = 0.042 and = 2 x 10 ^, and two different 
initial perturbation amplitude. Although at first the solution follows linear prediction. 



this trend reverses, consistently with Eq. (4.3 1, whose numerical solution is included in 



the figures. Such behavior seems to contradict linear theory, especially in Fig. 4.20 ^a), 
since the amplitude of perturbation at reversal clearly satisfies the linear ordering, 
1/ - /o|//o ~ 10"^ in this case. 

Subcritical solutions have also been identified with a different, two-species BB 
model, where thermal ion Landau d amping is included in a self-consistent way rather 
than as a constant input parameter [NLG^IO) . In this case, both resonant drive and 
resonant damping are present, and subcriticality is explained in terms of a nonlinear 
reduction of damping by particle trapping of the second population. An interesting 
feature is the existence of metastable, or subcritical steady-state solutions, which 
were not observed in single-species BB simulations. Since metastable modes exist for 
the two-species BB model in a regime consistent with BAEs in standard tokamak 
conditions [NGG"'"10] . subcritical BAEs can reasonably be thought as possible. It is 
then natural to raise a similar question, that is the existence of TAEs in a subcritical 
regime. Besides, an experimental device must take a path that goes through the 
linearly stable region as chirping emerges. Thus it is important to understand the 
mechanism of nonlinear drive in this regime. 

Here we show that, in our single-species BB model where external damping is 
treated as a fixed input parameter, subcriticality can be explained as an effect of 
holes and clumps. 
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Fig. 4.20: (a) and (c) Time-evolution of electric field amplitude obtained by linear 



theory (dotted line), by solving Eq. (4.3 ) (dashed line), and by solving the full-/ initial 
value problem (solid line), (b) and (d) Dynamic estimation of nonlinear growth rate 
(solid line). Dotted horizontal line is the value predicted by linear theory. 
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4.3.1 Mechanism 



Let us estimate the lowest-order correction to linear theory in the collisionless limit. 
Substituting E{x,t) = ujl{t) cos{ip)/k, where ip = kx — wt + <p(i), into the power 
balance Eq. (3.30), yields 



+ w / / cos(V')— dv + 7dWfc' = 0, 



(4.25) 



where the integral is taken over the resonant region. We decompose / in Fourier 
space, 



fix,v,t) = f{v,t) + [/„(i;,i)e"'^ + 



which simplifies the power balance, 

d^ ^ 
At 



ujTZ fidv, 



(4.26) 



(4.27) 



where TZ denotes the real part. Assuming the ordering / S> /i S> /2, the n = 1 
Fourier component of Vlasov equation yields 



dt 



where u = kv 



whose solution is 



Substituting the latter solution into Eq. (4.27) yields 



d^ 
dt 



LO ■ 



' ld)ojl 



2k 



n 



dif~fo) 
du 



(4.28) 



(4.29) 



dt', (4.30) 



which is actually an intermediary step in the derivation of Eq. (4.3 ). The latter equa- 



tion formally shows that small deviations in the velocity distribution, which are not 
taken into account in linear theory, could build up in time and eventually significantly 
modifies linear prediction. Physically, if we assume a fixed mode frequency, the drive 
is directly proportional to the velocity gradient in the resonant region. In Fig. [420l ;b) 
and (d), we estimate a nonlinear growth rate 7nl by measuring, at each time-step in 
the simulation, the slope of velocity distribution, averaged over a resonant region of 
width Aw = 4:UJii{t)/k. We also average this value in time with a window At — 50. 
We qualify 7nl as nonlinear, in the sense that it takes into account modifications of 
distribution function, although it is not estimated in a self-consistent way. At initial 
time, 7nl agrees with linear prediction, then it drops due to the flattening of the 
distribution. Afterwards, though there are cases where 7nl becomes positive, as in 
Fig. [420| ;b), there are also cases where it stays negative, as in Fig. [X20l ^d), which 
shows that we cannot assume a fixed mode frequency, since otherwise the wave would 
be stable in the latter case. 

Physically, subcritical instabilities can be explained by the following mechanism, 
which we illustrate in Fig. 



4.21 



by snapshots of / — /o in the neighborhood of resonant 
velocity. Even if the initial perturbation is small, some particles are trapped and create 
a seed island {jLt =10 — 20). As collisions try to recover the initial slope, small holes 
and clumps are created (7Lt — 30). As holes and clumps shift their central velocity 
{jLt — 150), they enhance free-energy extraction compared to a situation where the 
mode frequency would be fixed jBB95j . This interpretation is consistent with the fact 
that we only observe subcritical instabilities with frequency sweeping. We do observe 
subcritical solutions categorized as chaotic, but even chaotic solutions actually feature 
slight frequency sweeping (See Fig. |4.1[i)). 
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Fig. 4.21: Snapshots of deviation of velocity distribution from the initial one. Sim- 



ulation parameters are the same as for Fig. 4.20 |a) (we have increased to 8192 
in order to resolve small hole and clumps). Dotted line is a distribution such that 

ILO = Id- 



4.3.2 Initial amplitude threshold 

Ref. |BBC+99j gives a condition for subcritical instabilities to take off, as 



> 



1/2 

7lo 



(4.31) 



which is obtained by a dimensional analysis. We investigate this threshold with 
a series of Sf simulations with different initial perturbation amplitudes. A large 
value is taken for Ny, so that recurrence time is long enough (here Tr = 412/7lo)- 
In Fig. |4.22[ a), we observe a clear qualitative difference between solutions above 
and below uji,{t ~ 0)/jlo ~ 0-2. This threshold is consistent with numerical solu- 
tions of Eq. (4.3) shown in Fig. 4.22[ b). It is also in qualitative agreement with 
Eq. (4.311, which gives uji,{t = 0)/7lo ~ 0.15, with 7 — — O.Oll. The agreement be- 



tween Fig. 4.22 |a) and Fig. 4.22[ b) suggests that all information necessary to obtain 
an analytic prediction for the threshold, more precise than Eq. (4.31 ), are included in 
the reduced equation, Eq. 
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Fig. 4.22: Time-evolution of electric field amplitude in a subcritical regime, varying 
the initial perturbation amplitude, (a) (5/-C0BBLES simulations with 7^0 — 0.05, 
7rf = 0.062, Va = 10"-*, CFL = 0.65. (b) Numerical solutions of Eq. (lO). 
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Chapter 5 



Spectroscopic analysis of 
chirping TAEs 

As long as the background plasma parameters are not significantly changed, chirp- 
ing events in most tokamak experiments are quasi-periodic, with a quiescent phase 
between two chirping branches, which lasts a few milliseconds. It should be noted 
that this statement does not seem to apply to DIII-D |Hei95| . The spectrogram of a 
quasi-periodic chirping solution contains many pieces of information, which include 
chirping velocity, lifetime, quasi-period and asymmetry. This chapter shows how to 
use these information as diagnostics. We apply the BB model to AEs with frequency 
that sweeps only ~ 10% of the linear frequency, so that we can reasonably assume 
that no phase-space structure interacts with the bulk plasma. This situation is con- 
sistent with the Sf approach, where thermal populations are assumed adiabatic. In 
addition, the choice of Sf model removes complications of the full-/ approach that 
are due to effects of bulk. 

We consider TAEs in a time interval during which quasi-periodic, perturbative 
chirping is observed, and during which background plasma parameters are not signif- 
icantly changed. We assume a weak drive by energetic particles, such that chirping 
velocity is determined by a natural evolution of phase-space structures, rather than by 
a forcing from turbulent transport. We further assume that frequency shifting occurs 
well within the gap of Alfven continuum, so that chirping lifetime is determined by 
collision processes, rather than by continuum damping. 

In the spectrogram of magnetic field perturbations measured by a Mirnov coil at 
plasma edge, we extract the linear mode frequency /a, the average chirping veloc- 
ity, dSuP /dt, the maximum chirping lifetime, Tmax, and the average chirping period, 
Atchirp- The goal is to find the input parameters for 5f COBBLES simulations that 
are such that chirping shows similar features than those observed in experiments. 
When we compare simulation and experimental results, we simply re- normalize time 
by wa — '^t^Jpl, where /a is the TAE linear frequency. In this chapter, we adopt 
a convention of the EP literature, where growth rates and collision frequencies are 
given in percentage of wa- In the remaining of this work, all simulations are per- 
formed with Nx X 7V.U = 64 X 2048 grid points, and a time-step width Ai = 0.05, 
unless stated otherwise. A large number of grid points in velocity-space is necessary 
to avoid recurrence effect during quiescent phases between chirping events. 

5.1 Fitting procedure 

In the previous chapter, we saw that if we normalize time with the frequency of the 
mode, then chirping velocity, lifetime and period are dictated by the input parameters 
of the model, ^Lo^ Id, and, either j/q, either Vf and Vd- In the Krook case, we have 
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a 3-variables, 3-equations system, which we solve by a fitting procedure described in 
|5.1.1[ With drag and diffusion, there is one additional degree of freedom, hence the 
solution is not unique, but the boundaries of chirping regime limit the possible range 
of input parameters. 



Eq. (4.9) gives a relation between linear drive and external damping. 



2 1 dSuj^ 

7lo7. = (5.1) 

With the Krook model, chirping is limited to a range where 0.2 < 'fd/jw < 1-1 
|LIG09j . We found a similar constraint in our simulations with drag and diffusion, 
although a full scan of parameter space remains to be done. From this observation, 
in both cases, 7^0 is given within roughly 30% error, and 7^ within 50% error, by 



^Lo « 1-3 (t^^) , (5.2) 



1 dSuj^ \ ^ 
1 dSuj^ \ 9 



7^ « 0.7 (--^—) . (5.3) 



We refine these estimations in a manner that depends on the collision model we adopt. 
5.1.1 With Krook collisions 

The analysis described here aims at estimating the values of jloi Id, and Va for which 
the 6 J BB model fits experimental observations in terms of chirping characteristics. 



Eq. (4.22) yields the effective collision frequency, 

Va = —■ (5.4) 

^max 

Note that this effective collision frequency is meaningful only in the framework of a 
modelisation by the simple Krook operator of all dissipative processes, namely particle 
collisions, particle source and particle sink. Thus, Va can not be quantitatively com- 
pared with experimental measurements of collision frequency, unless particle source 



and sink terms are fully identified as well. Eqs. (5.1) and (5.4) form a system of 
two equations with three unknowns. The remaining unknown is found by fitting the 
chirping period. In our simulations, the chirping period is estimated by searching 
for the dominant frequency in the Fourier spectrum of electric field amplitude. To 
ensure a reasonable accuracy, simulations are performed for a time t 3> Atchirp- If 
the experiment belongs to a regime where /3 = 1, the above procedure is systematic. 
However, if /3 is significantly smaller than unity, an iterative procedure is needed, with 
a feedback between /3 and 7|q 7^. 



5.1.2 With drag and diffusion 



The analysis described here aims at estimating the values of 7^0, 7dj and Vd for 
which the 5f BB model fits experimental observations. Eqs. (4.24) and (5.1) form 
a system of two equations with four unknowns. The boundaries of chirping regime 
yield an estimation of Vd within 20% error. 



Vd 



1.2 



id 



1 dSuj' 



(5.5) 



On the one hand, it is shown in Ref. |LBS09| that for typical neutral beam-heated 
experiments, the ratio Vdjvf is of the order of unity. On the other hand, when 
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Fig. 5.1: Flowchart of chirping features fitting procedure in the drag/diffusion case. 
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> i^di we leave the periodic chirping regime for the dragging regime. Thus the 
relevant range for friction is Vf Vd- In this regime, Afchirp increases with decreasing 
Vf, 7, and increasing Vd- To obtain i^f and refine the above estimations, we need a 
two dimensional scan in (vf, Vd), where we search for solutions that fit the chirping 
period. In general, /3 7^ 1, and trial-and-errors are required to adjust chirping velocity 
to the experimental value. This procedure is summarized as a flowchart in Fig. |5.1[ 



5.2 Application to MAST 

The frequency sweeping mode in the MAST discharge #5568 between 64 and 73 ms 
has been identified as a global n = 1 TAB |PBG+04| . In the magnetic spectrogram 
shown in Fig. 5.2 'a), we measure /a = 111 kHz, d6uj'^/dt = 4.6 x 10~^, t„ 



' max 

0.35 X 10^, and Atchirp — 0.9 x 10^ (on average). 



5.2.1 With Krook collisions 



Eqs. (5.1) and (5.4) yield i^a — 3.1 x 10^^, and ^Lo^/ld = 1-5 x 10^^. However, the 
results of our analysis suggest that the plasma belongs to a regime where /3 = 0.85, 
hence we adjust the value of the product 7Lo^/7c^ to 1.5 x 10^^/0.85 = 1.8 x 10~^. 



A scan for this set of parameters is performed by changing the slope of the dis- 



tribution. Fig. 5.3 shows that the chirping quasi-period depends on in a roughly 
monotonous way. Note that the above scan needs to be performed on a relatively 
narrow range of 7^, since the limits of subcritical regime and non-chirping (chaotic) 
regime yield a first estimate as 7^ ^ 6 — 10% , 7^ ^ 2 — 7%. Here, the nonlinear 
stability threshold is defined as the largest value of 7^ for which the electric field 
amplitude vanishes in the time-asymptotic limit, independently of the initial pertur- 
bation amplitude; the chaotic regime is defined and categorized in a way described 



in Sec. 3.4.3 We observe that the two-points correlation of electric field amplitude 
decreases as the system approaches marginality. 

Fig. |5.2| [b) is the spectrogram for the simulation which is emphasized by a circle 
in Fig. |5.3[ The features of main chirping events agree with the experimental ob- 
servation. However, we observe a series of minor chirping events in between, which 
are absent in the experimental spectrogram. Another caveat is that only symmetric 
chirping is observed with the 5f BB model with Krook collisions and a linear velocity 
distribution. Thus, the application of this method with Krook collisions is restricted 
to symmetric or nearly-symmetric chirping experiments. Linear parameters estimated 
from this analysis are shown in Tab. |5.1[ Our analysis suggests that the TAE in this 
discharge is above marginal stability, with 7 ~ 7l. 



Eq. (4.8) yields a sa turated bo unce frequency as ~ 4.9%, which is satisfied in 
this simulation. In Ref. [PBG"'"04 , the peak amplitude of magnetic perturbation was 



determined with edge magnetic perturbation amplitudes and linear eigenfunctions 
calculated by the MISHKA code jMHKSQT] , and a scaling between Wf, and the mode 
amplitude was obtained with the HAGIS code. According to the latter analysis, 
Wf, ~ 5.4%, which is consistent with our result. However, this agreement is not 
enough to validate our procedure. 



5.2.2 With drag and diffusion 



We perform a first, rough scan in (vf, Vd) parameter space, assuming P — 1. Measur- 
ing average chirping velocity in repetitive chirping solutions yields an estimation of the 
correction parameter, /3 = 0.65. Then Eqs. (5.2), (5.3) and (5.5) yield 7^0 = 7.1*^%, 
7rf = 3.8^^%, and Vd = 1.6^° '^%. We perform a second, more careful scan, which con- 
sists of a series of 4 x 8 simulations in the domain (1.3% <Vd< 1.9%, 1 < Vd/vf < 8), 
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Fig. 5.2: (a) Spectrogram of magnetic fluctuations in MAST discharge #5568, ob- 
tained with a moving Fourier window of size 2 ms. (b) and (c) Spectrogram of the 
electric field, obtained with a moving Fourier window of size 0.6 ms, where the input 
parameters of the Sf BB model were chosen to fit the magnetic spectrogram for MAST 
discharge #5568. A solid curve shows the analytic prediction for chirping velocity and 
lifetime, (b) Krook collisions, correction parameter /3 — 0.85. (b) Friction-diffusion 
collisions, correction parameter f3 — 0.65. 
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Fig. 5.3: Scan of chirping quasi-period for the set of parameters corresponding to 
MAST discharge #5568. The hatched region correspond to a chaotic regime without 
frequency sweeping. The Hnear stabiUty threshold is indicated as 7 = 0, and the 
nonhnear stabihty threshold by a vertical dashed line. The chirping quasi-period 
agrees with the experiment for 7^ « 7.7%. Inset: zoom of the most relevant region. 
The spectrogram for the circled simulation is shown in Fig. |5.2|^b) . 
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Fig. 5.4: Time-evolution of the perturbation. The signal is filtered between 90 and 
130 kHz. The parameters of the simulation are shown in the second line of Tab. |5.l| 
For the simulation, to avoid hiding experimental data, we show the amplitude of 
perturbations (the envelope) instead of the perturbations themselves. Note the use 
of arbitrary units (we compare normalized quantities only) . 
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Tab. 5.1: Frequencies and growth rates estimated from the magnetic spectrogram of 
chirping TAEs, in percentage of the hnear frequency wa- 



Experiment MAST #5568 JT-60U E32359 

Coll. model Krook Friction-diffusion Krook Friction-diffusion 

"tZo 8T 9^4 9^8 

7l 7.7 10.5 8.5 9.2 

7rf 3.6 4.4 8.6 4.7 

Va 0.31 ... 0.25 

Vf ... 0.57 ... 0.36 

i/d ... 2.2 ... 1.7 

7 4.8 6.0 0.7 4.6 



where jlq and are constrained by Eqs. (4.24) and (5.1 ). The only repetitive chirp- 
ing solution with 800 < Atchirp < 1000 we found is shown in Fig. 5.2 'c). We verify that 
chirping features measured in this simulation fit the experiment, dSuj'^/dt = 4.7 x 10~^ 
(4.6 X 10-5 for up-chirping, 4.8 x 10"^ for down-chirping), r,nax = 0.35 x 10^ (0.36 x 10^ 
for up-chirping, 0.34 x 10^ for down-chirping), and A^chirp = 0.9 x 10"^. The estimated 
linear parameters are shown in Tab. |5.1[ Other solutions could be found with addi- 
tional simulations in the same domain, but the latter estimations are quite accurate 
because of the narrow range of periodic chirping regime. To validate this analysis, we 



compare the amplitude of perturbations in Fig. 5.4 Since the growth rate of chirping 



structure is neither 7 nor 7^, and the decay rate not simply 7^, but a function of sev- 
eral linear parameters, the agreement we obtain is not trivial (We measure a growth 
rate of 2.3%, and a decay rate of 0.6%). 

6.0%, 



Eq. (4. 



yields a saturated bounce frequency as ujb ~ 6.0%, which is also in 
agreement with the value of 5.4% estimated in Ref. [PBG+04) . 

In principle it should be possible to perform an independant estimation of collision 
frequencies using equilibrium parameters. However, some ambiguities in experimental 
data prevents such verification until more information are obtained. 



5.3 Application to JT-60U 



In JT-60U, TAEs are destabilized by a negative ion-based neutral beam (N-NB), 
which injects deuterons at Eb — 360 keV. In the discharge E32359, around t = 4.2 
s, fast-FS modes have been identified as m/n = 2/1 and 3/1 TAEs 'KKK+99'. In 



the spectrogram shown in Fig. 5.6 'a), we measure /a — 53 kHz, ddu^/dt — 6.3 x 
10"^, Tmax = 0.44 X lO'^, and Aichirp = 3 x 10^^ (on average). Frequency sweeps 
between 43 and 62 kHz. By comparing these frequencies with the SAW continuum 
gap in Fig |2.1[ it looks as if up-chirping was extending slightly into the continuum. 
However, the above continuum was estimated by assuming concentric circular flux 
surfaces, and a radial gradient of Shafranov shift can significantly increase the gap 
width. Unfortunately, current profiles are not available for this shot, forbidding more 
precise estimation of the continuous spectrum. However, there is a shot with similar 
background plasma parameters, E36378, for which current profiles are available and 
the continuum has been calculated |SKT+01| . In Fig. 5.5 we superpose the extent 



of chirping, measured in a magnetic spectrogram, to the SAW continuum calculated 
in the reference. At t = 3.65s, frequency sweeps between 53 and 61 kHz, while the 
continuum is broken between 41 and 76 kHz. At t = 3.75s, frequency sweeps between 
43 and 52 kHz, while the continuum is broken between 31 and 63 kHz. In both cases, 
the extent of chirping TAEs is well within the gap. Thus we can reasonably assume 
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Fig. 5.5: SAW continuum for n = 1, reproduced from Figs. 3 and 4 of Ref. [SKT+Ol] . 
normalized to ujq ~ Stt/q, to which we superpose a horizontal segment at the TAE 
linear frequency, and a vertical double-arrow, which represents the extent of chirping, 
(a) At t = 3.65, /o = 49kHz. (b) At t = 3.75, /q = 39kHz. 



that chirping lifetime is determined by collision processes, rather than by continuum 
damping. 



5.3.1 With Krook collisions 



Eqs. (5.11 and ( |5.4[ ) yield — 2.5 x lO^"^, and jLOy/ld — 1-8 x 10^ . However, a 
preliminary analysis suggests that the plasma belongs to a regime where /3 = 0.65, 
hence we adjust the value of the product 7^0^/7^ to 1.8 x 10^^/0.65 = 2.8 x 10^^. 

A scan for this set of parameters is given in Fig. |5.7[ The limits of subcritical 
regime and non-chirping regime yield 7l ^ 8 — 12%, 7(i 4 — 10%. Fig. |5.6[ b) is 
the spectrogram for the simulation which is emphasized by a circle in Fig. |5.7[ The 
features of main chirping events agree with the experimental observation. However, 
quiescent phases are replaced by many chirping events with amplitude slightly smaller 
than the major ones. Linear parameters estimated from this analysis are shown in 
Tab. |5.1| Our analysis suggests that the TAE in this discharge is marginally unsta- 



ble, with 7/7L ~ 0.1, even though 7^ < 7^;, which is not inconsistent with Eq. (3.59) 
since 7^0 > 7d- However, these values are inconsistent with estimations that take 
into account drag and diffusion processes. Since the following analysis with drag and 
diffusion shows much better agreement with the experiment, we imply that the Krook 
model is insufhcient to describe nonlinear features related to chirping repetition. 



5.3.2 With drag and diffusion 



A first scan in {vt^vj) yields an estimation of the correction parameter, /? = 0.85. 
Then Eqs. Q and ([5^ yield 7^0 = 10^^%, 7^ = 5±3%, and Vd = 1.7±°-3%. 

We perform a second scan, which consists of a series of 4 x 8 simulations in the domain 
(1.5% <Vd< 2.2%, 1 < Vd/vf < 8), where 71,0 and 7^ are constrained by Eqs. (4.24) 
and (5.1 ). The only repetitive chirping solution with 2500 < Atchirp < 3500 wc found 
c). We verify that chirping features measured in this simulation 
dSuj'^/dt = 7.1 X 10-^ (6.2 x IQ-^ for up-chirping, 7.9 x lO"'^ foj. 
Tmax = 0.45 X 10^ (0.47 X 10^ for up-chirping, 0.43 x 10^ for down- 



5.6 



is shown in Fig 
fit the experiment 
down-chirping) 
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Fig. 5.6: (a) Spectrogram of magnetic fluctuations during fast-FS modes in JT-60U 
discharge E32359, obtained with a moving Fourier window of size 2 ms. (b) and (c) 
Spectrogram of the electric field, where the input parameters of the 5f BB model were 
chosen to fit the magnetic spectrogram for JT-60U discharge E32359. A solid curve 
shows the analytic prediction for chirping velocity and lifetime, (b) Krook collisions, 
correction parameter /3 — 0.65. (c) Friction-diffusion collisions, correction parameter 
f3 = 0.85. 
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Fig. 5.7: Scan of the chirping quasi-period for the set of parameters corresponding 
to JT-60U discharge E32359. The region where 7l > 8.8%, where Afchhp < 2 x 
10'^, is not included in this plot. Both linear and nonlinear stability threshold are 



indicated. The chirping quasi-period agrees with the experiment for 
spectrogram for the circled simulation is shown in Fig. 5.6 'b). 



8.5%. The 




Fig. 5.8: Evolution of the perturbation during a single chirping event. The signal is 
filtered between 40 and 65 kHz. In these arbitrary units, 10^'' roughly corresponds to a 
noise level. The parameters of the simulation are shown in the fourth line of Tab. |5.l] 
For the simulation, to avoid hiding experimental data, we show the amplitude of 
perturbations (the enveloppe) instead of the perturbations themselves. 
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5.1 



chirping), and A^djirp = 3.1 x 10''. Estimated linear parameters are shown in Tab 
To vahdate this analysis, we compare the amplitude of perturbations in Fig. |5.8| We 
obtain a quantitative agreement, with a growth rate of 2.3%, and a decay rate of 
0.3%. 

As an independent validation, we now estimate the values of z/j and i^d from back- 
ground plasma parameters. In the discharge E32359 around t = 4.2 s, the resonant 
surface of the m/n = 2/1 and 3/1 TAE is located around r = 0.7 m. The magnetic 
shear is estimated from the q profile jGBC"'"Od] . 5* = 0.8. The deuteron plasma has 
the following characteristics, Bq ~ 1.2 T, Rq = 3.3 m, and the tangential radius of 
the N-NB is Rt = 2.6 m. At r = 0.7 m, Ue = 1.4 • lO'^^ m'^, and Tq = 0.75 
keV. We take into account carbon impurities with Z^ff — 2.7. With these equilibrium 
measurements, Eqs. ( 2.79||2.80 l yield Vfjujp^ ~ 1.2% and Vd/uJA = 1.7%. We obtain a 
quantitative agreement for Vd, which further validates our fitting procedure. However, 
Vf was underestimated by 70% compared to the latter estimation from background 
plasma parameters. Though error bars in experimental data may account for this 
discrepancy, it is also possible that our model misses some mechanism that would 
enhance friction. 

Note that electrons account for 99% of v'j^ which reflects a high Alfven velocity, 
while impurities account for 57% of v^, which is consistent with the fact that pitch- 
angle scattering is more effective with heavier particles. This observation suggests 
that impurities tend to reduce the ratio Vf/vdi which reduces the frequency of AE 
activity. 



5.4 Validity of fitting procedure 

Overall, our fitting procedure with drag and diffusion shows quantitative agreement 
for the time-evolution of perturbation, for an independent estimation of the linear 
drive in the case of MAST, and for an independent estimation of collision frequencies 
in the case of JT-60U. On the other hand, our fitting procedure with Krook colli- 
sions give similar linear rates only in MAST case, and in general we could not clearly 
reproduce periodic chirping with quiescent phases. Thus we infer that a fitting proce- 
dure of chirping features to estimate linear properties of background plasma is valid, 
provided that drag and diffusion are taken into account. 
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Chapter 6 

Conclusion 



Resonant interactions between NBI-induced far passing energetic ions and TAEs were 
described by a Hamiltonian perturbation in guiding-center angle-action variables. By 
reducing the latter Hamiltonian, a Fokker-Planck collision operator, and external 
damping processes to a new, 2D phase-space, we showed how the input parameters 
of the BE model are quantitatively related to plasma parameters and characteristics 
of the TAE. 

We developed a kinetic code to solve the initial-value BB model. COBBLES was 
verified by checking conservation properties, benchmarked against numerical simula- 
tions in former works, and validated against well-known linear and nonlinear theories. 
The feasibility of long-time simulations for a I0W-7L distribution hinged upon quick 
convergence and strong numerical stability of the CIP-CSL scheme. 

With both Sf and full-/ versions of COBBLES, we explored theory in several 
parameter regimes, namely steady-state, periodic, chaotic, chirping, and subcritical. 
We performed a scan of the nonlinear evolution of electric field in the whole (7^, 
Va) parameter space. A new diagnostics allowed ns to identify the chirping region. 
Although holes and clumps were not expected to appear when < 0.4 7/,, we found 
that the frequency sweeping region can expand to a low external dissipation regime, 
around ~ 0.2 7^. 

Numerical results display good quantitative agreements with theory in several 
regimes. The limits of validity correspond to the assumptions used in theory. A 
perturbative numerical approach which doesn't take into account the kinetic response 
of the bulk may feature spurious agreement outside of the validity limits when the 
resonant region reaches a non negligible portion of bulk particles. 

We found a regime of quasi-periodic chirping with both Krook and drag/diffusion 
collision operators. Since quantitative agreement with theory suggests the predictabil- 
ity of nonlinear chirping characteristics based on fundamental linear kinetic param- 
eters, the latter may be estimated in the opposite way from chirping data in exper- 
iments. More precisely, chirping velocity and lifetime yield two relations among 7^, 
7ci, and collision frequencies; and a fitting of Atchirp yields an estimation of remaining 
unknowns. Note that major advantages of this technique are 1. kinetic parameters 
in the core of the plasma estimated only from a spectrogram of magnetic fluctua- 
tions measured at the edge, without expensive kinetic/MHD calculations nor detailed 
core diagnostics, and 2. unified treatments of supercritical and subcritical AEs. We 
showed that diffusion is essential to reproduce quiescent phases observed in experi- 
ments between chirping events, and drag is required to observe repetitive chirping, 
since drag allows the distribution to recover in a reasonable time. We confronted this 
procedure by analyzing AEs on MAST and JT-60U. We found quantitative agreement 
with measured magnetic fluctuations for the growth and decay of chirping structures, 
quantitative agreement with an estimation of the linear drive based on measured 
magnetic fluctuation amplitude in the case of MAST, and qualitative agreement with 
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collision frequencies estimated from experimental background measurements in the 
case of JT-60U. In the latter estimations, impurities, which were not included in es- 
timations of Ref. |LBS09| . account for the main part of velocity diffusion. In both 
cases, our analysis suggests that TAEs belong to a regime which is relatively far from 
marginal stability, where total growth rate is of the same order as linear drive. 

Accurate estimations of linear drive, damping rate, and overall growth rate, open 
the way to promising advanced numerical and experimental investigations. 

1. Applying the same procedure to dozens of available chirping TAE data would 
produce a database from which useful trends could be extracted. 

2. The estimated growth rates and collision frequencies could be translated in 
terms of spatial gradient of energetic-particle distribution, assuming a slowing- 
down distribution in energy, and in terms of coefficients of the 3D Fokker-Planck 
collision operator. This work would yield all input parameters of drift-kinetic 
perturbative 3D simulations that assume a fixed and arbitrary external damp- 
ing, such as in the HAGIS code, where drag and diffusion is being implemented. 

3. Since 7 is roughly proportional to ujdw f + ndp^f, where W is the energy, it 
should be possible to control the growth rate by varying the total number of 
injected high-energy ions. A possible scenario is to start from a NBI experiment 
with chirping TAEs, and reproduce the same conditions except for different NBI 
power Pnbi at the time of chirping. With a scan in Pnbi, dependency of 7 on 
Pnbi and existence of subcritical TAEs could be investigated. 
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Appendix A 



Validity limit of Lie transform 
perturbation theory 



A.l Test problem 

To illustrate the improvements of Lie transform perturbation theory compared to 
classical perturbation theory, let us apply both formulations to a simple mathematical 
problem, 

^1 = -I - (A.2) 

where e is a small parameter, with initial conditions x|^^q = and y|^^o ~ These 
are the equations of motion of an Hamiltonian system, with Hamiltonian 

9 9 

H^iV) "y + y + 4^'' = ho{x,y) + €hi{x). (A. 3) 

In action-angle variables z = ((^, J), which are defined by 



X = V2JsinC, (A.4) 

y = V2JcosC, (A.5) 

the unperturbed motion is described by 

m - i, (A.6) 

J{t) = K, (A.7) 

and the perturbed system by the total Hamiltonian, 

hiC, J) = J + eJ^ sin-*(C). (A.8) 

A.2 Classical perturbation theory 

Classical perturbation theory yields a n*'^-order solution of Hamilton's equations. 



dzo dzi dzn f I \ dh 



dt ' ^ dt dt \ -1 I dz 



, (A.9) 

z=Zo + ...+e" 1 z„_i 



after an expansion in the small parameter e, z = ZQ + ezi +e Z2 + . . .. Fig. A.l (a) and 



(c) show the evolution of action and total energy given by the application Eq. (A.9) 
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Fig. A.l: Evolution of action and energy, (a) and (c), classic perturbation theory, (b) 
and (d), Lie transform theory. In both cases, e = 0.1. 



at successive orders. These results are compared with the exact solution obtained by 
solving the original Hamilton's equations with a fourth-order Rungc-Kutta scheme 
[PTVF92] with a very small time-step width. At small times, the solution given by 
this method is close to the exact solution, but we rapidly observe a strong secularity 
effect that makes the average action spuriously shrink or expand. 



A. 3 Lie perturbation theory 



Let us now treat the same problem, this time with Lie perturbation theory. The 
fundamental one- form, 7 = 70 + £71 + e^72 + . • ., can be separated into its equilibrium 
part and the perturbation. 



70 = J rfC - ho{J) dt, 

71 = -hi{C,J)dt, 



(A.IO) 
(A.ll) 



and 7„ = for n > 2. 

To simplify the fundamental one-form, we change variables from z to Z = (^, /). 
Following Lie perturbation theory, the zeroth-order one-form and the time coordinate 
remain unchanged, Tq = 70, and T = t. Applying Eqs. (2.37), (2.39), and (2.38), we 



find the new one-form, gauge scalar, and generating functions, which yield a solution 
which is valid up to n*'^-order in e. We choose initial conditions such that orbits have 
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a same energy level. At first order, 



1 + 



/ = K. 



(A.12) 
(A.13) 



This solution can be obtained in terms of the original coordinates by using the pull- 
back operator, 



c 

J 



4 

/2 r 



2sin(2e) - ^sin(40 



^ cos(4C) - 2 cos(2^) 



At second order, 

i = 

I = 

and the pull-back operator is, 
/ 



4 64 



64 256 



c 



4 

512 



2sin(20 - |Sin(4e) 



(A.14) 
(A.15) 

(A.16) 
(A.17) 



e-— [sin(80 - 16sin(6C) - 4sin(40 -I- 400sin(2^)] -I- o(e^)(A.18) 



J = I 



64 



cos(4^) - 2cos(2^) 



[2cos(60 + 6cos(4C) - 42cos(20 + 17] + 0(6^). (A.19) 



Fig. A.l (b) and (d) show action and total energy for analytic solutions, Eq. (A.12 
A.15 1 and ( A.16|[ATT9 ), with e = 0.1. There is good agreement with the exact solution 
without any secularity effect. 

Fig. |A.2| shows the relative error in total energy conservation against the small 
parameter e. It is confirmed that the error is one order higher than the order of 
perturbation analysis. 
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Fig. A. 2: Relative error in total energy conservation for solutions obtained with Lie 
perturbation theory at t = 12.57r. For small e, solid, dotted, and dashed lines scale 
as e, e^, and e'^, respectively. 
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Appendix B 



Validity limit of gyrokinetics 
in slab geometry 

Most of linear and nonlinear numerical investigations of TAEs are based on a hy- 
brid model where energetic particles follow gyrokinetic equations, while background 
thermal populations are described by a MHD model |TS98| . Such simulations can 
exhibit chirping that is qualitatively similar to experiments |TSTI03] . Ultimately, the 
results we obtain in this thesis with our ID approach should be compared with global 
gyrokinetic simulations - actually, this should have been an intermediary step before 
carrying out comparisons with experiments. 

This appendix summarizes some results of gyrokinetic theory, which could be used 
in future works. Note that the following results are not used within the core of this 
manuscript. The goal is to investigate a break-down of gyrokinetics when time-scale 
and length-scale of the perturbation approach the gyromotion time-scale and length- 
scale, and when the amplitude of perturbation approaches the equilibrium. Since we 
used Lie-transform perturbation theory to derive the set of gyrokinetic equations, this 
provides a second illustration of this technique. 

B.l Review of gyrokinetics 

Let us consider the evolution of a distribution / in 6D (3D position, 3D velocity) 
phase-space z, given by a Vlasov equation, 

dtf - {h, /}, = 0, (B.l) 

where h is the Hamiltonian, and {}z arc Poisson brackets. This description can be 
simplified by removing the irrelevant fast gyromotion from the complete motion of 
charged particles, thus reducing dimensionality to 5D. There are two ways to proceed, 

• averaging Vlasov equation over a gyroperiod, 

dt if)^ - {{h, = 0, (B.2) 

where ^ is the cyclotronic angle, or 

• finding a coordinate transformation from z to Z, where / and h are changed to 
F and H such that 

dtF - {H, F}^ = 0, (B.3) 

and H does not depend on ^. Lie transform provides a modern way of deriving 
gyrokinetic equations in the latter fashion, which preserves the Hamiltonian structure 
of the initial problem. These conservation properties are essential for robust long- 
time computations. Fig. |B.1| is a schematic summary of the modern derivation of 
gyrokinetic equations. 
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Fig. B.l: Schematic summary of the modern derivation of gyrokinetic equations for 
an electrostatic perturbation. 
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B.2 Single particle orbit in slab geometry 



We develop a numerical simulation to solve the equations of motion for a single ion 
in both particle coordinates, guiding-center coordinates, and gyrocenter coordinates. 
The aim is to compare 3 ways of computing the orbit of a single particle. 

We consider an electrostatic perturbation ip of frequency and wave vector w and 
fc, respectively. For the gyrokinetic treatment to be valid, the gyrokinetic ordering 
|Hah88j . 

« 1, (B.4) 

i i^i '^J. -^i 

k±p, « 1, (B.5) 

has to be satisfied. Here is the ion cyclotronic frequency, Ti the ion temperature, 
pi the ion Larmor radius, and L„ a scale-length of variation of equilibrium plasma. 

As a test problem, we assume an homogeneous magnetic field Bo = Bob and a 
time-independant, single-wave electrostatic perturbation, 

(p{x) — ifQCOs{k-x). (B.6) 
Then, the gyrokinetic ordering is reduced to ^ ~ ^ ^ I, and k±pi w 1. 

B.2.1 Equations of motion 

We normalize length, time, velocity, and electric potential, to Debye length Xjj, inverse 
ion plasma frequency ion thermal velocity vxi, and the ratio T,;/e respectively. 

In particle position coordinates 

The equations of motion expressed in particle position coordinates are given by New- 
ton's equations with Lorentz's force, 

X = V, (B.7) 
■i) = E + n^v X b. (B.8) 

In guiding-center coordinates 

We recall the guiding-center variables Z = {R, U, M, ^), 

R = X e^a, (B.9) 
U = v-b, (B.IO) 
M = (B.ll) 



^ EE tan"' . (B.12) 

\v-e2j 

The equations of motion expressed in guiding-center coordinates |Lit83| are given by 

it = v\\b + n-^bxVip, (B.13) 

if = -b-Vip, (B.14) 

M = -n,,d^^, (B.15) 

i = n,,{l + dM^)- (B.16) 
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In gyrocenter coordinates 



In Ref . jHah88] , a complete set of gyrokinetic equations is derived in the homogeneous 

magnetic field case, based on Lie-transform perturbation theory. In the gyrocenter 
variables Z = {R, U, M,^), the equations of motion are 

il = Ub + n-^bxV if) , (B.17) 

TJ = -b-V{>fi), (B.18) 

11 = 0, (B.19) 

t = + (B.20) 



The first equation yields the velocity of the gyrocenter as a sum of parallel velocity 
and E X B velocity. The second one shows the effect on parallel acceleration of the 
electric field. The third one confirms that the new magnetic momentum is an exact 
invariant. The last one can be ignored in practical gyrokinetic simulations. However, 
it is evolved in our simulation since it is needed to retrieve guiding-center coordinates 
from gyrocenter coordinates and compare physical quantities expressed in different 
sets of coordinates. 



B.2.2 Numerical results 
Test problem 

Let the plasma be such that ujpi = 10 GHz, vtl = 10^ m.s^^ and Xd = 10 /im, with 
a magnetic field amplitude B = IT. Then fici = 96 MHz and pi = 1.04 mm, or 
flci = 9.6 X 10"'^ and pi = 104 in normalized units. We define two small parameters 
EE = and Ek = k\\/k±. In our test problem, we choose ee = 0.1, k± = 1/Pi, 
fc|| = O.lfcj^. The initial conditions are f = 0, and v± = = 1, so that M{t = 0) = 
0.5. 

This set of parameters satisfies the gyrokinetic ordering, 

CB-efe « 1, (B.21) 
k±p, ~ 1. (B.22) 



Comparison of numerical results in several coordinate systems 

To compare numerical results between themselves, we transform each coordinate sys- 
tem into guiding-center coordinates. For the particle coordinates we define 

Rxv = x-^a{C) (B.23) 

^ 'ci 

f/,, = vb (B.24) 

= ^ (B.25) 

C.. = r, (B.26) 

where 



Vj_ 



^/vei^ + ve2'^, (B.27) 
C = tan-^ {vei,ve2) . (B.28) 

Guiding-center coordinates are taken as the reference, Zgc = Z. 
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Fig. B.2: Time-evolution of a particle position in an uniform magnetic field with an 
electrostatic perturbation. Here, guiding-center trajectory overlaps the trajectories of 
both pull-back of gyrocenter, and push-forward of particle position. Simulation with 
e = 0.1 and At ^ 10'^ T^. 



Finally the guiding-center coordinates are related with the gyrocenter coordinates 
Z by the pull-back transform, 



^gy 


— V$ X b 

- R eE 2. ' 

Ct 


(B.29) 


^gy 


- 


(B.30) 


Mgy 




(B.31) 


^gy 




(B.32) 



where we have explicitly written the small parameters as reminders of the order of 
each term. Note that the generating vector for the parallel velocity is one order higher 
than the others. 

In order to compare results obtained by the gyrokinetic code with the true orbit of 
the particle and its guiding-center, we firstly run a numerical simulation with a very 
small time-step. At = lO^'^ Tci, where Td = 27T/ilci is a gyration period. Fig. B.2 
shows the orbit of the particle calculated in position/ velocity coordinates, and the 
orbit of the guiding-center, which overlaps with orbits of both the push-forward of 
the particle position and the pull-back of the gyrocenter. Then we choose a practical 



time-step width, and compare, in Fig. B.4 the guiding center Hamiltonian calculated 
in the three coordinates system. 



H. 



H. 



gy 



M. 



1 



gc -r 2^gc + <^E fiRgc, Mgc, ^gcj , 
1, 



+ ^U'^y + '^E 'PiRgy , Mgy , ^gy) . 



(B.33) 
(B.34) 
(B.35) 
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Fig. B.3: (a) Magnetic momentum, (b) Gyroangle perturbation (difference between ^ 
and the unperturbed solution — ^dt)- Simulation with e = 0.1 and At = lO^'^ Td. 
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Fig. B.4: Error in energy conservation. The time-step width of these simulations is 
At = Te,/2. 
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Fig. B.5: Error in energy conservation aX t — 100 T^i, against time-step width. A 
solid line follows a law, and dashed lines follow a At^ law. Two points marked 
[A) show a break-up of gyrokinetics close to cyclotronic period. 



We observe a secularity effect in the numerical error of particle position. The gyro- 
center Hamiltonian, 



(B.36) 



is trivially perfectly conserved. 
Analysis of the error 

Firstly, we analyse the numerical error. Our code uses a fourth-order Runge-Kutta 
scheme jPTVF92] . As a consequence, the error between numerical and analytic solu- 
tions is proportionnal to At^. Fig. B.5 shows the variation of the relative error in the 
conservation of guiding-center and gyrocenter Hamiltonians, calculated after a fixed, 
long time (hundred gyration periods), with increasing time-step. This plot is consis- 
tent with the prediction of numerical error for this numerical scheme. Note that the 
error for each coordinate is proportional to At^ , save the magnetic momentum calcu- 
lated in gyrocenter coordinates, which is trivially exactly conserved for any time-step 
width. The curve of gyrokinetic-calculated error crosses the other ones at a fraction 
of the cyclotronic period, meaning that a gyrokinetic computation becomes rapidly 
more efficient. However, the points marked (A) show that even when the gyrocenter 
Hamiltonian is precisely conserved, the conservation of guiding-center Hamiltonian 
breaks-down when the time-step width is close to the cyclotronic period. At ~ T^. 
This is caused by a lack of accuracy of the pull-back transformation of coordinates, 
coming from a lack of accuracy for large time-step width in the computation of 

Fig. |B.6| shows the variation of the relative error in the conservation of guiding- 
center and gyrocenter Hamiltonians with the small parameters ce and e^, for a very 
small time-step width. Fig. |B.7| shows the variation of the relative error in the conser- 
vation of guiding-center and gyrocenter Hamiltonian with the small parameter e, when 
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k\\/k±. A solid line follows ej. law. Simulation parameters are At — O.OlTd, and 
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£k = = £■ To predict the variation of the relative error in the conservation of the 
Hamiltonian (B.35), let us calculate its difference with the gyrocenter Hamiltonian, 
whose error in conservation does not depend on the small parameters, 



H 



Mgy - M + - (U'gy - [/') + eE viZgy) - CE (ifi) (Z), (B.37) 



-EE (p{Z) - - {Ugy + U) EEek-^ + eE (p{Z) 




<ci 




The latter result is consistent with both Fig. |B.6| and Fig. |B.7| 

B.3 Conclusion 

Using a full set of gyrokinetic equations, which is valid up to first order in the small 
parameter e, in the case of an homogeneous magnetic field and a small electrostatic 
perturbation that satisfies the gyrokinetic ordering, we computed a single charged- 
particle orbit. We confirmed that the error due to the gyrokinetic reduction is of 
the order of . This simple test problem shows the accuracy of gyrokinetic equations 
compared to a direct computation of the particle orbit in position/velocity or guiding- 
center coordinates. We note that in practical gyrokinetic codes, information about 
the gyroangle is discarded, because of its inaccuracy and its irrelevancy. 
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Appendix C 

The CIP-CSL scheme in 2D 
phase-space 

To solve the general ID advection equation, 

dtF + udxF = 0, (C.l) 

in a conservative form, we evolve both F and its integrated value R. -F" (i?") is the 
continuous interpolated solution for F {R) at a time t„ = nAt. F^ is the discrete 
value of F" at the grid point of coordinate Xm = mAX, and 

R^^ = /'"''F«(A)dA. (C.2) 



F" is an array which contains the values F^ for every point m on a grid point of 
coordinate A„. We define the ID algorithm CIPCSLIdIu, F", i?", R"+\ A, 

AA, At) as follows. 
For each m, 



• Define Ap„ = Xm — uAt 



Find the grid point km satisfying A^^ < Ap„ < A^^+i, and define (Am) 



Define the coefficients 



Pk + Fk +1 2i?" 
AA2 AA3 



(C.3) 



• Advect R, 
where 

• Advect F, 



+ FP , , 3R2 , , 

Vk^ = + (C.4) 



= (Am)' + (Am)' + Fll (Xm) ; (C.6) 



3</)fc^ (Am)' + (Am) + i^fcl- (C.7) 

This ID algorithm is extended to the 2D phase-space {x,v) in the following way. 
is the value of the distribution / at the grid point of coordinates Xi = iAx, 
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= v^ia + j Aw, at a time t„ = nAt. We define the density within a cell and the 
line densities, 

P{x,v)Ax<lv, (C.8) 



J Xi ^ 

r{x,v,)dx, (C.9) 
r{x,,v)Av. (C.IO) 

j 

The first advection in the a;-direction is performed by calling successively, for each j, 

CIPCSLlD{vj, a*, p", p*, a:. Ax, At/2), 
CIPCSLlD{v,, r, /*, trj, aj, X, Ax, At/2), 

with tJj = {vj + Vj+i)/2. Similarly, the advection in the i;-direction is performed by 
calling successively, for each i, 

CIPCSLlD{qEjm, a*, a**, p* , p** , v, Av, At), 
CIPCSLlD{qEjm, f*, f**, a*, a**, v, Av, At), 

with Ei = (Ei + Ei^i)/2. Then we repeat the advection in the a;-direction, 

CIPCSLlDivj, a**, p**, X, Ax, At/2), 

CIPCSLlD{vj, /**, /"+!, a**, crj+i, x, Ax, At/2). 

Here, subscripts * and were used to designate intermediate steps between i„ and 
tn+i- To avoid spurious leakage of particles, we impose a zero flux at the velocity 
boundaries by setting 

(A,=i) = (A,.^J = 0. (C.ll) 



Boundary conditions in the velocity distribution are illustrated in Fig. C.l^a), 
which shows the discretization on Ny = 8 grid points of a Maxwellian distribution 
f{v), and its interpolation polynomial F. Notice two buffer points j — and j = 
Ny + 1. Fig. |C.l[ b) shows the integral function, 

d{v) = /' f{v')dv', (C.12) 



and its interpolation polynomial D. Discretization of p, which is defined as 

p{v,) = f{v')dv', (C.13) 



is included. Here we chose a small number of grid points and small cut-off velocities 
to emphasize border effects. For realistic parameters, the discrepancies between / 
and F, and between d and D are negligible. 
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Fig. C.l: Initial Maxwellian distribution discretized on — 8 grid points, (a) 
Particle distribution /, interpolation polynomial F, and discretized values fj. (b) 
Density distribution d, interpolation polynomial D, and discretized values of grid- 
point density pj. Inset: zoom on the first points. 
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Appendix D 



Nonlinear regime 
categorization algorithm 



The algorithm we developed to categorize the behavior of each numerical solution 
is an improved version of the algorithm designed by Vann [VDR^OS) . with an other 
way of sorting out chaotic from periodic behavior, a new diagnostics to distinguish 
chirping from merely chaotic solutions, and where we take into account numerical 
issues that appear when 7^ is smaller. 

The categorization is based on an analysis of the time-series of normalized electric 
field energy, A{t) = £{t)/£o, with £0 = / 2. Firstly, we define the global minimum 
as ^gm = iiiiii{^(^)}o<t<tn,ax: where tmax is the time-duration of the simulation. Then 
we drop the initial transient phase, and extract a time window < t < tmax over 
which A{t) is sampled at a rate Atg. Here, imin is an estimation of the time-duration 
of the transient phase. Within this window, we define the mean value (A), maximum 
Ajnax = ma.x{A{t)} and minimum ^min = min{v4(t)}, the oscillation amplitude AA = 
^max — ^min, SLild the local minima (maxima) as points where A{t) is smaller (larger) 
than A{t — Ats) and A{t + Atg). As a measure of the periodicity, we compute the 
two-point correlation function as 

(i,(r,r')i,+i(r,T')) 



^^^^ = ^T. . ^ ..I.. 

"^.=0 (A,(r,r')2) (A,+i(r,r')2) 



where, for each correlation window size r, 

i,(T, r') = A(W - J T - t') - (A(<max - 3 ^)) , (D-2) 

'Ti- = (^max ~ ^min)/''' is the number of periods included inside the total time-window, 
and angular brackets represent a time-average over a period, 

Air.r')) ^\ ^'^^ A{T,T')dT' . (D.3) 

The overall correlation i?o is defined as the maximum of -R(t) for r > Tc, where Tc 
is the shortest period such that R{t^ < 0. In other words, Rq is the normalized 
amplitude of the peak in the two-point correlation function corresponding to the 
dominant frequency. This measure, which is used to sort out chaotic from periodic 
behavior, is different from the one given in Ref. |VDR+03| . We also provide a criterion 
to sort out chirping solutions. We compute the frequency power spectrum of the 
time-evolution of the electric field at some position, E{x = 0,<), and normalize the 
amplitude of the spectrum to its maximum value. At each time, we extract the largest 
and smallest frequency for which the amplitude in the power spectrum is significant. 
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i.e. larger than a threshold e^. Then, we define Aw^ax as the maximum difference 
between the largest and smallest frequency, normalized to the plasma frequency. We 
then proceed to the following decision tree: 

1. IF Agn, < eo THEN damped 

2. ELSE IF Awmax > £7 THEN chirping 

3. ELSE IF (A) < ei THEN damped 

4. ELSE IF A{t) is monotonic AND AA/ (A) < 62 THEN steady-state 

5. ELSE IF A{t) is monotonically decreasing (zero local extrema) THEN damped 

6. ELSE IF each minima (maxima) is larger (smaller) than the former OR A A/ (A) < 
£3 OR AA < £4 THEN steady-state 

7. ELSE IF i?o > 1 - £5 THEN periodic 

8. ELSE IF the number of extrema is not less than four THEN chaotic, 

where are thresholds that must be carefully adjusted. Special care is taken in 
adjusting empirically so that frequency splitting is not mistaken for frequency 
sweeping. 

In this decision tree, the logical test [T] is an addition to the decision tree given 
in Ref. [VDR~'~03j . For damped solutions, as the electric field becomes small, the 
particles experience free streaming, leading to spurious recurrence effects after half a 
recurrence time = n/kAv. Free streaming occurs when EqTr < Av, and cq 

is chosen to reflect this condition. For the benchmark in |3.4.3[ this logical step is 
switched off as the recurrence effect is less problematic for shorter-time simulations. 
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Appendix E 



Chirping features analysis 
algorithm 



Numerical validation of chirping theory, and development of semi-empirical laws for 
chirping features, both hinge upon systematic quantification of chirping velocity, life- 
time, and period or quasi-period. 

We developed an algorithm to measure features of major chirping events in both 
full-/ and Sf COBBLES simulations, for both Krook and drag/diffusion collisions, 
except in the dragging regime. The starting point is a spectrogram P{ti,ujj), such as 



plotted in Fig. 5.6 ^b) or (c), which is obtained from time-series of the electric field at 
a; = with a moving Fourier window of size At, and the linear frequency uj, which 
is obtained from one of our linear analysis tools described in Sec. |3.2[ Here, i is an 
index of the median time of a Fourier window, and j is an index of the frequency, 
where i^j+i — = 27r/At. Let us consider upward chirping only. The first step is 
to identify major chirping branches. A simple approach would be to extract from 
the spectrogram, at each time step, the local maximum, and try to fit a square root 
function following maxima whose amplitude exceeds some threshold. However, two or 
more chirping branches including minor events can easily be confused as one (Minor 
and major events are not well separated categories, it is difficult to set a threshold to 
isolate major events). Trying to start from the tip of a branch and follow it down by 
jumping to the closest maximum fails for the same reason. The problem with these 
procedures is that, when we extract maxima, we loose useful information, namely the 
amplitude of these maxima. The following approach makes use of some additional 
information. It requires a rough estimation of chirping velocity ~ {dSuP' / dtY^"^ ^ 
and a rough estimation of chirping lifetime r^ax; both of which are easily obtained 
from theory using input parameters of the simulation. First, we initialize a test 
chirping velocity V to . For each time-step ti„ of the spectrogram, the algorithm 
takes the sum Si^^iV) of P{ti,u}j), where i starts from iq, and is such that ojj{ti) 
roughly follows a square-root law with a coefficient given by V . We repeat this 
procedure for each iq, and for test chirping velocities O.IF" <V < \QV'^ . A maximum 
of Sig {V) should indicate the initial time ti^ and velocity V of a major chirping event. 
We found this is the case most of the time for data that are analyzed in this work. 
Simultaneously, we extract the maximum Pmax reached by P during this chirping 
event. The second step is to identify the tip of chirping. Using the estimated chirping 
velocity and r^^xi deduce its approximate location i^jp and jjjp. The algorithm, 
starting from well above j^jp and decrementing j, searches around i^jp for the first 
value of P above a threshold which is given by e^^Pmax- This should correspond to 
the tip of the considered chirping event, itip and jtip- Since the algorithm is sometimes 
mistaken, we go on decrementing j, while following the first local maximum to the 
left (left meaning smaller time here). When we reach a j such that LOj k, oj, we 
check that i is close to Then the lifetime is given by ti^^ — ti^. Otherwise the 
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chirping event is discarded. Repeating the above steps while avoiding the time region 
ti^ < t < ti^^^ yields velocity and lifetime of a significant number of major chirping 
events. Finally the chirping period is simply obtained by extracting the maximum in 
the Fourier spectrum of electric field time-series. The whole procedure requires only a 
few seconds on one typical laptop processor, even for long-time simulations {t ~ 10^). 

The same algorithm is used to measure chirping features in experimental data. In 
this case, the starting point is a spectrogram of magnetic perturbations measured by 



a Mirnov coil at the edge of the plasma, such as Fig. 5.6 [a,), the linear frequency oj and 



a rough estimation of chirping velocity, and of chirping lifetime, which are directly 
measured on the spectrogram. 
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